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Abstract 

We present a new general relativistic code for hydrodynamical supernova simulations with neutrino 
transport in spherical and azimuthal symmetry (ID and 2D, respectively). The code is a combination 
of the CoCoNuT hydro module, which is a Riemann-solver based, high-resolution shock-capturing 
method, and the three-flavor, fully energy dependent VERTEX scheme for the transport of massless 
neutrinos. VERTEX integrates the coupled neutrino energy and momentum equations with a variable 
\ Eddington factor closure computed from a model Boltzmann equation and uses the "ray-by-ray plus" 

, approximation in 2D, assuming the neutrino distribution to be axially symmetric around the radial 

direction at every point in space, and thus the neutrino flux to be radial. Our spacetime treatment 
I employs the Arnowitt-Deser-Misncr (ADM) 3-1-1 formalism with the conformal flatness condition 

(CFC) for the spatial three-metric. This approach is exact for the ID case and has previously been 
shown to yield very accurate results for spherical and rotational stellar core collapse. We introduce 
new formulations of the energy equation to improve total energy conservation in relativistic and 
, Newtonian hydro simulations with grid-based Eulerian finite-volume codes. Moreover, a modified 

version of the VERTEX scheme is developed that simultaneously conserves energy and lepton number 
in the neutrino transport with better accuracy and higher numerical stability in the high-energy tail 
of the spectrum. To verify our code, we conduct a series of tests in spherical symmetry, including 
O ■ a detailed comparison with published results of the collapse, shock formation, shock breakout, and 

^ ' accretion phases. Long-time simulations of proto-neutron star cooling until several seconds after 

c/3 I core bounce both demonstrate the robustness of the new CoCoNuT-VERTEX code and show the 

^ ■ approximate treatment of relativistic effects by means of an effective relativistic gravitational potential 

as in PROMETHEUS- VERTEX to be remarkably accurate in spherical symmetry. 
. Subject headings: supernovae: general — neutrinos — radiative transfer — hydrodynamics — relativity — 

^ ' methods: numerical 
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1. INTRODUCTION 



The collapse and explosion of massive stars has been an active topic in astrophysical research for decades, and has 
provided a major impetus for advances in computational relativistic astrophysics, in particular ever since the first 
numerical simulations of core collapse were conducted in the 1960s (Colgate & White 1966). The general outline of 
this rich astrophysical phenomenon is fairly well established by now: After the collapse and bounce of the central iron 
core to a compact neutron star {M/R 0.1 . . .0.2 in relativistic units) a shock wave is formed, which stalls a few 
ms after bounce. Afterwards, accretion onto the newly-formed neutron star continues typically for a few hundreds 
of ms until the shock is re-energized either by neutrino heating (the "delayed neutrino-driven mechanism" , Wilson 
1985), probably aided by multi-dimensional hydrodynamical instabilities like convection Epstein 1979; Herant, Benz, 
^ & Colgate 1992; Herant et al. 1994; Burrows, Hayes, & Fryxeh 1995; Janka & Muller 1995, 1996; Buras et al. 
■ - ' 2006b, 2006a; Bruenn et al. 2006; Burrows et al. 2006a and the standing accretion-shock instability (SASI; Blondin, 
Mezzacappa, & DeMarino 2003), or by some alternative mechanism, such as MHD effects (e.g. Burrows et al., 2007 and 
references therein), or energy deposition by acoustic waves (Burrows et al. 2006b). The ongoing quest for the explosion 
mechanism, with various research groups pursuing a number of different avenues, is a testimony to the complexity 
of the problem, which involves not only multi-dimensional (magneto-)hydrodynamics, but also kinetic theory (i.e. 
neutrino transport), nuclear and particle physics, and general relativistic (GR) gravity. Even for the few hundreds 
of ms from the collapse to the onset of the explosion, there is no unified approach which fully accounts for all these 
intricacies. Instead, a number of quite disparate approaches has been used in the past to investigate different aspects 
of core collapse supernovae, all of which are arguably deficient in the sense that they sacrifice some essential physical 
ingredients in favour of others. 

Thus, there have been repeated attempts to solve the problem of the explosion mechanism by means of an accurate 
treatment of neutrino transport in spherical symmetry, combined with reasonably accurate microphysics (i.e. the 
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nuclear equation of state, nuclear burning, and neutrino-matter interactions). These endeavours culminated in the 
development of solvers for the neutrino Boltzmann equation (Mezzacappa & Messer 1999; Burrows ct al. 2000a; 
Liebendorfer at al. 2001b, 2004; Rampp & Janka 2002; Yamada, Janka, & Suzuki 1999); but the exclusion of multi- 
dimensional hydrodynamical effects is obviously a serious shortcoming. On the other hand, the restriction to spherical 
symmetry allows for a full GR treatment with relatively little additional effort: Relativistic supernova simulations with 
grey or multi-group flux-limited diffusion have been conducted since the 1980s in ID (Baron et al. 1989; Bruenn, De 
Nisco, & Mezzacappa 2001), and relativistic Boltzmann transport became feasible after another few years (Yamada 
ct al. 1999; Liebendorfer et al. 2001b, 2004). As expected. GR effects were shown to be anything but negligible: 
Comparisons with the Newtonian case (see, e.g. Bruenn et al. 2001) revealed a sizable influence of GR on the mass 
of the inner core at shock formation, the initial shock strength, the long-term evolution of the shock position, the 
compactness of the proto-neutron star, and the neutrino luminosities and energies. 

Another complementary approach has focused primarily on the interplay between neutrino heating/cooling and multi- 
dimensional (magneto-)hydrodynamic effects, which has turned out to be no less important for the understanding of the 
explosion mechanism than the; proper treatment of neutrino transport and general relativity: Convective instabilities 
were found to operate both inside the proto-neutron star (Epstein 1979) and the "hot-bubble region" behind the stalled 
shock (Bethe 1990). Hot-bubble convection in particular was shown to be helpful for increasing the neutrino heating 
efficiency already in the 1990s (Herant et al. 1992, 1994; Burrows et al. 1995; Janka & Miiller 1995, 1996), and another 
recently identified hydrodynamic instability of the shock, known as the "standing accretion shock instability" (SASI; 
Blondin et al. 2003), was found to have a similar effect. Nowadays, the most advanced computational tools used in 
this approach combine multi-dimensional Newtonian hydrodynamics with an up-to-date treatment of the microphysics 
with multi-group neutrino transport methods of similar sophistication as in the ID case (Livne et al. 2004; Buras et al. 
2006b; Bruenn et al. 2006; Liebendorfer, Whitehouse, & Fischer 2009). But some weak spots remain to be criticized: 
First, multi-dimensional neutrino transport has not been treated in its full complexity yet; the approaches now in use 
either assume the neutrino flux vector to point in the radial direction (the ray-by-ray approach; Buras ct al. 2006b; 
Bruenn ct al. 2006), neglect energy redistribution in phase space (Livne et al. 2004; Liebendorfer et al. 2009), or 
involve even more radical approximations. Second, all serioiis attempts to model the development of multi-dimensional 
instabilities during the post-bounce phase, while treating the neutrino transport with reasonable accuracy, have been 
carried out in the framework of Newtonian hydrodynamics, which is rather unsatisfactory considering the importance 
of GR effects established in ID studies. So far, attempts to include relativistic corrections have been limited to 
the use of modified gravitational potentials that can be heuristically derived from the Tolman-Oppenheimer-Volkov 
equation of stellar structure (Rampp & Janka 2002; Marek et al. 2006). Although this method can be justified on 
the basis of comparisons with results from relativistic ID neutrino transport codes and multi-dimensional relativistic 
hydrodynamics codes without neutrino transport (Liebendorfer et al. 2005; Marek et al. 2006; Miiller, Dimmelmeier, 
& Miiller 2008), it comes along with several drawbacks: As the equations of Newtonian hydrodynamics are left 
unaltered, it may lead to incorrect predictions for some dynamical properties of the neutron star such as the oscillation 
eigenfrequencies (Miiller et al. 2008), which may also result in a distortion of the gravitational wave spectrum. 
Moreover, the "effective potential approach" ignores the effects of relativistic kinematics and still allows superluminal 
velocities - which have actually occurred in Newtonian simulations of the explosions of 0-Ne-Mg cores. Finally, it 
should be noted that, different from the purely Newtonian case, a conservation law for the total (i.e. kinetic, internal, 
and potential) energy cannot be formulated when a modified gravitational potential is used. As a consequence, it is 
difficult to exclude the possibility of non-physical energy loss or generation, in particular as cumulative effect over long 
evolution periods. 

On the other hand, the last few years have seen the development of genuinely relativistic hydrodynamics codes 
for multi-dimensional simulations of core collapse, often based on a rather crude treatment of the microphysics, and 
without any up-to-date scheme for the neutrino transport. Typically, simple analytic equations of state (EoS), such 
as the hybrid EoS of Janka, Zwerger, & Monchmeyer (1993) were employed until recent years; and the polytropic 
models for the progenitor had only a limited resemblance to iron cores from stellar evolution calculations. These 
limitations have naturally led the numerical relativity community to focus on the one hand on the only phases of the 
supernova for which such severe approximations may be deemed sufficient, i.e. the collapse and bounce of the core, 
and on the other hand on those issues such as black hole formation (Sekiguchi & Shibata 2005), gravitational wave 
emission from rotational core collapse (Dimmelmeier, Font, & Miiller 2002a, 2002b; Shibata & Sekiguchi 2004), and 
the development of non-axisymmetric instabilities (Shibata & Sekiguchi 2005) that do not lend themselves readily to 
a Newtonian analysis, at least not with the same accuracy. Obviously, such an approach, which neglects physics that 
is unquestionably of paramount importance for the stellar core collapse problem, is beset with large imcertainties. 
However, the first steps towards a more sophisticated treatment of the microphysics in core collapse supernovae were 
taken by Dimmelmeier et al. (2007) and Ott et al. (2007a), who studied the collapse of rotating iron cores with a 
nuclear equation of state and a simple parameterized scheme for the deleptonization during collapse (Liebendorfer 
2005). The fact that their findings are rather different from earlier studies with more radical approximations - they 
exclusively find gravitational wave signals of the so-called type I (Zwerger & Miiller 1997) - underlines the need for a 
reliable treatment of the microphysics also in the field of relativistic stellar core collapse. 

Clearly, the most satisfactory approach to modeling supernovae would be to sacrifice as little of the essential ingre- 
dients (neutrino transport, equation of state and other microphysics, general relativity, multi-dimensional hydrody- 
namics) as possible in the numerical models. In this paper, we therefore describe a relativistic generalization of the 
Newtonian ray-by-ray-plus method of Buras et al. (2006b) for the neutrino transport, which we haved combined with 
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the so-called conformal flatness condition (CFG) or Isenberg- Wilson-Mathews approximation (Isenberg 1978; Wilson, 
Mathews, & Marronetti 1996) for the gravitational field equations. The method presented here has been implemented 
by building on the existing radiation transport code VERTEX (Rampp & Janka 2002; Buras et al. 2006b) and the 
relativistic hydrodynamics code CoCoNuT (Dimmelmeier et al. 2002a, 2005), and has already been applied to 2D 
simulations of the collapse and post-bounce phase, which will be the subject of future papers. At this stage, however, 
we shall confine ourselves to presenting thorough tests of our new VERTEX-CoCoNuT code in spherical symmetry. 

This paper is organized as follows: In Sec. [5J we discuss the governing equations of hydrodynamics and neutrino 
transport ( "neutrino hydrodynamics" ) for the ray-by-ray-plus approximation, and comment on the their implementa- 
tion in CoCoNuT and VERTEX. Tests of the new code in spherical symmetry are presented in Sec.[3l most notably 
a detailed comparison of ID simulations of the 15 Mq progenitor model of Woosley & Weaver (1995) with VERTEX- 
CoCoNuT, the relativistic Boltzmann code AGILE of Liebendorfer et al. (2004), and the Newtonian VERTEX code 
of Rampp & Janka (2002), in which GR effects are taken into account approximately by using a modified gravitational 
potential. In this context, the issue of energy and lepton number conservation also receives particular attention. In 
order to demonstrate the ability of VERTEX-GoGoNuT to evolve supernova models stably well beyond 1 s after core 
bounce, we also present results for a ID explosion model of an 0-Ne-Mg core (Nomoto 1984, 1987; Kitaura, Janka, & 
Hillebrandt 2006). 

Further technical details of our numerical implementation, which may also be of relevance for similar radiation 
hydrodynamics codes, are given in the Appendices: In App. |A] we address the problem of total energy conservation 
for a self-gravitating fluid in Newtonian and relativistic hydrodynamics, and in App. [B]we outline an efficient finite- 
difference scheme for the neutrino moment equations that guarantees the simultaneous conservation of energy and 
lepton number. 

Note that geometrized units are used throughout Sec. [21 i.e. both the speed of light and the gravitational constant 
are set to unity: G = c = 1. Greek indices run from to 3, Latin indices from 1 to 3. 

2. GENERAL RELATIVISTIC NEUTRINO RADIATION HYDRODYNAMICS 

2.1. Hydrodynamics 

In general relativity, the hydrodynamic evolution of a perfect fiuid is governed by two conservation equations for the 
baryonic rest mass current J** and the stress-energy tensor T^^^ , 

= 0, V.T^'^ = 0, (1) 

where V i, denotes the covariant derivative. For a perfect fiuid, and T^" can be expressed in terms of the baryon 
rest- mass density p in the local fiuid frame, the four- velocity w^, the pressure P, the specific enthalpy h = l + € + P/p 
(where e is the specific internal energy), and the four- metric g^^, 

= pu^\ Tf"" = phut^u" + Pg^"'. (2) 

For the metric, we adopt the Arnowitt-Deser-Misner (ADM) 3 + 1 formalism (Lichnerovicz 1944) of general relativity 
to foliate the spacetime into spacelike hypersurfaces. In this approach, the four-dimensional line element ds^ — 
g^i^ dx'^ dx'^ is written as follows, 

ds2 ^ ^a^di^ + (da;* + P'dt) [dx^ + (3^dt) , (3) 

where a is the lapse function, /3* is the shift vector, and 7^ is the induced three-metric on each hypersurface. Using 
this decomposition of the four-metric, the equations of hydrodynamics can be formulated in the frame of an Eulerian 
observer'^ moving orthogonally to the spacelike hypersurfaces, i.e. with a four-velocity = (a~^, q;~^/3*). As in 
Banyuls et al. (1997), we introduce the following conserved variables^, 

D = pW, = phW^v\ T^phW^-P-D, (4) 

^ Since the terminology in the existing literature on relativistic hydrodynamics and radiative transfer is not uniform, we add some 
remarks on this issue here. The "frame" of an observer (or a class of observers) is an orthonormal tetrad field whose time-like base vector 
coincides with the four-velocity of the observer. We shall use the terms "Eulerian frame" (as in Banyuls et al. 1997) or "lab frame" (as in 
Mihalas & Weibel Mihalas 1984) interchangeably to denote the frame of an Eulerian observer as described in the main body of the text. 
It is important to note that such an Eulerian observer is not a fixed-coordinate observer unless the shift vector vanishes. 

D, Si and t can be obtained form and T^" using and the projection operator ±tJ onto the three-hypersurfaces: D = n^J^^ , 
= ±l,n^Ti^'', and r = n^n^Tf'^ - D. 
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where = /(avP) + /a is the three-velocity in the Eulerian frame, and W — 1/ \/l — ViV^ is the corresponding 
Lorentz factor. The equations of GR hydrodynamics in flux-conservative form then read, 



dt dx^ 2^ " dxi V dt 

~9 



c 




dxf' ^ ^^"'J ' \ dt ) 



^^Q^lna ^ ypypo 



c 

(7) 

Here, is defined as = — /3*/q;, and g and 7 are the determinants of the four-metric and the three-metric 
"fij, respectively. Different from the purely hydrodynamic case discussed by Banyuls et al. (1997), the right-hand sides 
(RHSs) of Eqs. (ini and ([7]) contain source terms for the exchange of momentum and energy with neutrinos (which can 
be related to angular moments of the collision integral and are therefore denoted by a subscript "C") in addition to 
the gravitational and geometric source terms. These source terms will be specified in Sec. 12.3.51 

Eqs. ([5]) to ([7]) must be supplemented by an additional equation expressing the conservation of electron number in 
the absence of weak interactions (which will be discussed separately) , 



dt ^ dx' ~ \ dt J 



(8) 

c 

where Ye is the electron fraction (number of electrons minus number of positrons per baryon). The change of the 
electron number due to neutrino emission or absorption is accounted for by the source term on the RHS, in a similar 
fashion as energy and momentum exchange with neutrinos in Eqs. © and ([7]). 

Since we are dealing with a multi-component fluid whose composition is not uniquely determined by the thermo- 
dynamic state (p, r. Ye) unless nuclear statistical equilibrium (NSE) applies, additional conservation equations for the 
mass fractions of protons, neutrons, a-particles and heavier nuclei are also needed, 

dt ^ dx- ' 

Finally, the pressure and the composition (in the NSE regime) in Eqs. (|6I7I9P have to be provided by an equation of 
state (EoS). 

We rely on the CoCoNuT code (Dimmelmeier et al. 2002a, 2005) to solve the equations of relativistic hydrodynamics 
in spherical polar coordinates by means of a high-resolution shock-capturing (HRSC) scheme, employing piecewise 
parabolic (PPM) reconstruction (Colella & Woodward 1984), an approximate Riemann solver, and second-order Runge- 
Kutta timestepping. Our method of choice for the Riemann solver is a new hybrid HLLC/HLLE scheme in CoCoNuT 
along the lines of Quirk (1994), combining the high resolution of the HLLC solver (Toro, Spruce, & Speares 1994; 
Mignone & Bodo 2005) with the robustness of the HLLE solver (Harten, Lax, & van Leer 1983; Einfeldt 1988) against 
odd-even decoupling near grid-aligned shocks (Quirk 1994; Lion 2000; Kifonidis et al. 2003; Sutherland, Bisset, & 
Bicknell 2003). In order to treat the advection equations for the nuclear mass fractions, we have added a simplified 
version of the consistent multi-fluid advection (CMA) method by Plewa & Miiller (1999). Moreover, the computational 
performance of the code has been improved significantly by enhancing parallelism in the axisymmetric (2D) mode. 

Most notably, however, our most recent simulations also use a reformulated energy equation instead of Eq. ([7|) to 
improve total energy conservation. The new scheme, which is described in detail in App. \^ is designed to alleviate a 
well-known problem with the gravitational source term in the energy equation; While there is a global conservation 
law for the ADM energy (or the total energy in the Newtonian case), the gravitational source term in the standard 
form of the energy equation Q used in Dimmelmeier et al. (2002a, 2005, 2008) does not have the form of a flux 
divergence, and therefore a standard finite-volume discretization will lead to a numerical violation of global energy 
conservation. Nonetheless, as we demonstrate in App. \^ the gravitational source term can partly be absorbed in a 
flux divergence term to reduce this conservation error considerably. In the case of a stationary spacetime, or in the 
Newtonian limit, it is even possible to achieve total energy conservation up to machine precision. 



2.2. Metric Equations - The Conformal Flatness Condition 

The CoCoNuT code solves Einstein's field equation using the conformal flatness condition (CFG), which was 
introduced by Isenberg (1978) and first used in a pseudo-dynamical context by Wilson et al. (1996). In the GFG 
approximation, the spatial three-metric is assumed to be conformally flat^, i.e. it is obtained from the flat three- 
metric^ 7ij by a conformal transformation 7^- = 0*7ij, where <p is the conformal factor, thus reducing the number 
of independent metric quantities to five. As a consequence, the four constraint equations for the metric in the ADM 

^ We note in passing that this assumption can always be fulfilled by an appropriate gauge choice fo r a s pherically symmetric spacetime, 
and no approximations to the ADM equations are actually made by using CFC in that case. See Sec. 12.41 for more details. 
® Thus, in spherical polar coordinates we have 7ij = diag (l, , sin^ 9) . 
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formalism, combined with a slicing condition, determine the metric completely. In our case, we use maximal slicing, 
i. e. we require the trace of the extrinsic curvature K'^^ to vanish: K = Kl = 0. We then obtain a set of five non-linear 
elliptic equations for a, 0, and /?% 

A<^--2^<^^(^£;+^^g^), (10) 
A(a0)==-2™0^ S + 25+^-^^ , (11) 

\ lOTT / 

kp' = levra^^S^ + 2(j)^°K''\/, (^^^ - ^V^V,/3^ (12) 

where A and V are the Laplace and covariant derivative operators for a flat three-space. The total matter-energy 
density E = t + D (for r and D see Sec. p.ip ). the three-momentum density Si and the trace S = ^ijS^^ of the spatial 
components Sij of the stress-energy tensor (as measured by an Eulerian observer) appear in the non-linear source 
terms on the RHS, and comprise contributions both from matter and (neutrino) radiation. For the matter component 
we obtain 

E = phW^-P, (13) 
S = phW^v^+5P, (14) 
S'^phW^v\ (15) 

while for neutrinos, 

E = S = 4TrW^ { J + 2vrH + vfK) , (16) 

Si='iTT(t)"^W^ [H (l + w2) + Vr (J + K)] , (17) 
^2 = ^3 = 0, (18) 

in our approximation. Here, J, H, and K denote the zeroth, first and second moments of the neutrino radiation 
intensity in the comoving orthonormal frame (see Sec. 12.31) . and u,- — (jp'v^ is the radial velocity component in the 
orthonormal basis of an Eulerian observer. 

We solve the CFC equations ((TO)) . (ITlT) . and (fT2)) in the reformulated version of Cordero- Carrion et al. (2009) using 
a fixed point iteration scheme as described in Dimmelmeier et al. (2005); i.e. during each iteration step, the non-linear 
source terms are updated, and the Laplace operators on the left-hand side (LHS) of the equations are inverted by means 
of a multipole expansion (cf. Miiller & Steinmetz 1995). While this method suffers from a somewhat low convergence 
rate compared to the spectral nonlinear Poisson solver of Dimmelmeier et al. (2005), the underlying algorithm is more 
easily amenable to parallelization, which is a critical issue for multi-dimensional neutrino transport simulations. 

2.3. Neutrino Transport in CFC Spacetime 
2.3.1. Overview 

Up to now, all self-consistent numerical simulations of neutrino transport in core collapse supernovae are based on 
the assumptions i) that a semi-classical treatment of kinetic theory is applicable, and ii) that neutrino oscillations 
need not be taken into account (see Cardall, 2008; Yamada, 2000; Strack & Burrows, 2005 for a discussion of the 
more general case). In our approach to the transport problem, we retain those assumptions, although they may be 
invalid under certain specific circumstances (cf. Sec. 12. 4p . Under this proviso, the evolution of the invariant neutrino 
distribution function / (the number of neutrinos per phase-space volume d'^x d'^p) is governed by the relativistic 
Boltzmann equation for massless or ultra-relativistic particles, which reads (see e.g. Lindquist, 1966; Ehlers, 1971; 
Stewart, 1971) aj- j tj. a/- 

D^— + — — = £ fl9) 

in any coordinate basis. Here, / and the collision term £ depend on the spacetime coordinates and the four- 
momentum vector pf^ (obeying the mass-shell constraint p^^p^ = 0) as measured in the associated holonomic base 9^; 
A is the affine parameter. 

In general, Eq. ()19p describes a six-dimensional time-evolution problem, which can be scaled down considerably by 
using the ray-by-ray-plus approximation (Buras et al. 2006b) in the following manner: First, the equations of neutrino 
transport are formulated for the spherically symmetric case fSec. 12.3.'^ to 12. 3. 5t . thus reducing the dimensionality of 
the individual transport problems to be solved from six to three. By working with a finite number (i.e. two) of angular 
moments of the Boltzmann equation, yet another dimension can be eliminated. The additional space dimensions are 
then taken into account by solving these equations independently on different radial "rays" (corresponding to the 
angular zones of the polar grid) with direction unit vector n, assuming rotational symmetry around n for the neutrino 
distribution function. In addition, certain terms from the full transport equations in 2D (or 3D in general) are taken 
into account to avoid unphysical behaviour in the optically thick regime (see Buras et al., 2006b and Sec. 12.3.5)) . We 
emphasize that the ray-by-ray-plus approximation does not mean that the transport of neutrinos is considered to be 
only in the radial direction. Instead it is assumed that the neutrino distribution function is axially symmetric around 
the radial direction, and that the neutrino pressure tensor is diagonal. 
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2.3.2. Boltzmann Equation in Spherical Symmetry 

The precise form of Eq. (jl9p in a spherically symmetric spacetime has been discussed in various places in the literature 
(Lindquist 1966; Castor 1972; Schinder 1988; Mezzacappa & Matzner 1989); it depends on the adopted gauge and 
slicing conditions, and on the tetrad basis used in momentum space. In our case, the spacetime metric is of the CFC 
type with = and p"^ = ^ 0, 

ds^ = -a(t,rfAt^ + (j}{t,rf \{Ar + p'{t,r) dtf+r^dO'^+r'^sm^edifA . (20) 

Furthermore, we carry out a transformation to the orthonornial frame conioving with the fluid, in which the collision 
integral £ can be expressed most conveniently; the expressions for the collision integral provided in App. A of Rampp 
& Janka (2002) thus remain valid. With the velocity field chosen as 



(</)-2i,^(t,r),0,0). 



(21) 



a possible choice for the basis four- vectors of such a frame in terms of the coordinate derivatives 9t, dr, dg and (see 
Hawking & Ellis, 1973; Misner, Thorne, & Wheeler, 1973; Straumann, 2004 for this notation) is given by 



ei = 

63 = 



a-^VrWdt + W \ 



a~^f3^Vr) dr, 



(22) 
(23) 
(24) 
(25) 



In spherical symmetry, the distribution function / depends only on t, r, the comoving frame energy e = p • gq (where 
p is the neutrino four- momentum) and the angle cosine = p • ei/p • Gq, and obeys the following PDE: 



W 



9/ 
dt 

da 
^ dr 
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dt 



dr 
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dvr 



dt 
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dr 
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dVr 
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dr 
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dr 



(26) 



where v = [i 



^ = 1 + [iVr (bearing in mind that c = 1) and Tp — 1 ~ fi^. It should be noted that our 



result is analytically equivalent to the one obtained by Mezzacappa & Matzner (1989) for the same gauge and slicing 
conditions. The transfer equation for the radiation intensity I = h~^c~^£^ f can be obtained by exploiting the relation 
e^df /de = d {fs^) /de ~ 3e^/, and is omitted here for the sake of brevity. 

2.3.3. Exact Moment Equations for CFC Metric m Spherical Symmetry 

Multiplying Eq. ([26]) by h~^c~^e^ and taking the zeroth and first angular moments yields the moment equations in 
a spherically symmetric CFC spacetime, resulting in 



dW{J + VrH 



dt 
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for the momentum equation. Here, the angular moments J, H , K and L of the specific intensity are given by 
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(28) 



(29) 



while C^^-' and C*^^^ are the corresponding zeroth and first moment of the collision integral. Note that Eqs. ([Tf)) and 
are written in terms of the densitized moments J, H, K and L (where X = for any quantity X) which 

contain the square root of the determinant of the three-metric ^/^ = (jfir^ svaO as an additional factor, in order to 
better reflect the underlying law of energy conservation. The virtue of such a formulation is more apparent in the 
equations governing the evolution of the neutrino number density and flux. 
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where {J^, H, JC, C,C} = {e^^J, e^^H, e^^K, e^^L, e^^C}. Apart from the colhsion integral, no source terms appear in 
the neutrino number equation pop. Since the evolved quantity W i^J + VrhC) is proportional to the spectral neutrino 
number density A'oui in the Eulerian (coordinate/lab) frame, 

47r 

N,^^=—W{J^Vrn), (32) 
c 

a finite-volume discretization of Eq. pop exactly conserves the total neutrino number measured in the laboratory 
(coordinate) frame, f f An 



j j de /77Veui = ^ Jd^x de (J 



VrU) , 



(33) 



in the absence of neutrino-matter interactions. As a consequence, numerical conservation of lepton number can be 
achieved with an appropriate conservative treatment of the fiux terms in the neutrino moment equations and the 
source term for the electron fraction. 

With a few minor exceptions, the numerical solution of the neutrino moment equations in our code proceeds exactly 
as described by Rampp & Janka (2002). Equations (P7|) and (|28l) are written in a conservative finite-difference form 
using second-order space discretization on a staggered mesh and backward differencing in time; only the neutrino 



advection terms 



d_ 

dr 



dr 



a 



(34) 



are discretized with upwind differences in space and central differences in time. The higher moments K and L are 
computed from the variable Eddington factors Jk and /l obtained from the solution of a model Boltzmann equation 
(cf. Sec. 12.3.41) as K — fxJ and L = f^J where necessary. 

Our numerical scheme for the moment equations differs from that of Rampp & Janka (2002) in one important 
respect. Rampp & Janka (2002) also solve Eqs. (1501) and (pij) for the neutrino number density and flux in addition 
to Eqs. (j27p and (I28p in order to guarantee both lepton number conservation and energy conservation. However, this 
procedure has two important drawbacks: Solving another pair of equations is computationally expensive, and, more 
importantly, the numerical solution can become inconsistent in the sense that { J,H) = {eJ^^sT-L) no longer holds. In 
the high-energy tail of the spectrum Ji/mathcalJi and Hi /Hi in a given zone in energy space can even assume values 
outside the zone boundaries. In order to circumvent this problem, we have developed a special second-order finite- 
difference representation for the terms governing advection in energy space accounting for gravitational redshift and 
Doppler shift. This allows us to guarantee both energy and lepton number conservation without solving the equations 
for neutrino number density and flux. Since our improved treatment of the redshift and Doppler terms is rather 
involved, we give a complete description of the algorithm in App. [BJ 

Finally, the moment equations are supplemented by two additional equations for the operator-split update of the 
specific internal energy and the electron fraction of the matter due to neutrino-matter interactions. 
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where the factor aW ^ takes care of the conversion from proper time to coordinate time. The resulting non-linear 
system of equations for J, H, e and Ye is then solved by Newton- Raphson iteration. As in Rampp & Janka (2002), 
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we treat fi and r neutrinos and their antiparticles as a single species, and solve the corresponding moment equations 
separately from the (fully coupled) moment equations for the transport of electron neutrinos and antineutrinos. 

2.3.4. Model BoHzmann Equation 

In the variable Eddington factor technique, the closure relations required for Eqs. (|27p and (|28|) are extracted from 
the formal solution of a model Boltzmann equation, with the moments of the neutrino distribution function occurring 
in the collision integral taken from the solution of the moment equations in an iterative procedure (Mihalas & Weibel 
Mihalas 1984; Burrows et al. 2000b; Rampp & Janka 2002). Since the model Boltzmann equation is only used 
to compute normalized moments of the radiation intensity, it is sufficient to consider a strongly simplified equation 
version of Eq. (|26p. We therefore neglect terms containing the shift vector or derivatives of the metric functions a and 
(j) in Eq. ([26| . as well as higher-order terms 0{v'^) in the radial velocity. Furthermore, we set = + Vr ^ fi and 
^ = 1 + iJ,Vr — )• 1 in the terms for advection in energy space. However, the radial advection term is kept exactly as in 
the full equation. The model transfer equation for the radiation intensity thus obtained reads 
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dfj. 
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(37) 



Here, C [I] denotes the collision term in the transfer equation, which is related to the collision term £ in the Boltzmann 
equation ([25]) by C [I] — e'^h^^c^'^(i[f]. Aside from the appearance of prefactors containing the metric functions a 
and (j) and of the radial shift (3^ in the advection term, Eq. (1571) is identical to the flat-space transfer equation in the 
0(w/c)-approximation. Therefore, the same solution method as in Rampp & Janka (2002) can be employed, i.e. the 
radial advection term and the energy derivatives in Eq. p7p are treated by a conservative interpolation procedure and 
a time-explicit upwind scheme respectively, while the tangent-ray method (Yorke 1980; Mihalas & Weibel Mihalas 
1984) is used for the remaining terms. 



2.3.5. Source Terms 

Once the solution of the moment equations is known, the source terms Qy^ Qe, and Qm for electron number, energy, 
and momentum in the local fluid frame can be computed from the zcroth and flrst angular moments C^*^^ and C^^^ of 
the collision integral. 



Qy,= 
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(38) 



(39) 



Qm = -4:Tt / de 



E 



(40) 



To obtain the correct transformation to the Eulerian frame (in which the equations of hydrodynamics are solved), 
we consider the equations of electron number conservation and energy-momentum conservation in their covariant 
form. Source terms due to neutrino interactions must enter as scalar and four-vector terms on the right-side in these 
equations: 

ipY.u'^) = s, V.T'^'^ = (41) 

Here, s and are the scalar electron number density source term and the covariant energy-momentum vector source 
term. Since the transformation of scalars and four- vectors is trivial, we only need to express s and q^^ in terms of Qy^ , 
Qe, and Qm in the comoving frame to obtain the correct source terms in any desired frame: 

s = Qy, , = QeBo + QmBi . (42) 

After carrying out a Lorentz transformation^ back to the Eulerian frame and recasting the resulting equations into 
conservative form, we obtain the following source terms for the conserved Eulerian quantities that appear in Eqs. ([6]), 

Since we consider a purely radial velocity field at this stage, the Lorentz transformation of a four-vector from the comoving frame to 
the Eulerian frame simply reads: (q'^, q^, q^, q^) — > W{q^ + Vrq^, q^ + Vrq^, q^, q^). 
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and dll): 



V dt 



VjW Me + Qm) , 



- VrQu) , 



(43) 



(44) 



(45) 



Note that we have multiphed the source terms by the metric ^/j, because it is the densitized variables ^J^Si^ J^t and 
^J^DY(, that are actuahy evolved in Eqs. ([5]), © and ([5]). This is possible because the metric factor is obviously 
not changed by neutrino interactions. 

2.3.6. Ray-by-ray-plus Approach for Axisymmetric Problems 

As noted by Rampp & Janka (2002), a treatment of neutrino transport in spherical symmetry (as described in 
Secs. l2.3^ to l2.3.5)) can be extended to the multi-dimensional case by assuming that all the involved physical quantities 
depend only parametrically on the coordinates 6 and </> (using spherical polar coordinates), while neglecting any lateral 
derivatives in the full three-dimensional moment equations. However, Buras et al. (2006b) pointed out that the lateral 
advection of neutrinos in the optically thick regime and the non-radial components of the neutrino pressure gradient 
need to be taken into account to avoid unphysical convective activity in the proto-neutron star. In general relativity 
the additional advection terms in the moment equations are. 



f d^W {J + VrH) 

V dt 
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09 
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(46) 
(47) 



which are solved in an operator-split step before the radial transport sweep. We also include terms for the compressional 
heating of the neutrinos due to lateral motions, which are added to the moment equations (j27p and (j28[) when 
performing the radial transport sweep. 
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(49) 



It should be noted that gravitational redshift associated with lateral motion is still neglected here; however, this is a 
tiny effect in typical supernova environments due to the relatively small asphericity of the gravitational field. 

The lateral component of the neutrino pressure gradient appears as a source term for the component 5*2 of the fluid 
momentum density, 

( d^S2 \ ^dP, 



\ dt 



lat 



de 



(50) 



where we assume that the neutrino pressure is given by Pi, = 4/37r J in the optically thick regime. Since the acceleration 
four-vector due to the neutrino pressure gradient must be orthogonal to the four- velocity of the fluid, the source term 
in the energy equation (work exerted on the fluid) is given by 



( d^T \ 

V dt A,, 
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,dP 
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(51) 



Finally, contrary to the Newtonian case, the covariant source vector has to be transformed to the Eulerian frame 
by a Lorentz boost which includes the non-radial velocity components. Qe and Qm now enter differently into the 
source term for 5*1, and additional source terms for S2 and 6*3 appear. 
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(53) 
(54) 
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While the new terms in Qm arc of order 0(v'^) and ean usually be neglected, all the terms in Qe are of order 0{v) 
and need to be included to ensure the correct evolution of the specific internal energy e. This is particularly crucial 
for rapidly rotating configurations, where the rotational velocity can reach a significant fraction of c. 

2.4. Quality of our Approximation 

Our relativistic ray-by-ray-plus method entails several approximations, which fall into three categories, viz., approx- 
imations in the treatment i) of the microphysics (EoS, neutrino interactions, etc.), ii) of the neutrino transport, and 
iii) of the gravitational field equations. Concerning the microphysics, we rely on the same up-to-date description as 
the non-relativistic VERTEX-PROMETHEUS code (Rampp & Janka 2002; Buras et al. 2006b), which includes a 
number of additional neutrino reactions compared to the "standard" set of neutrino opacities used by Bruenn (1985) 
(see Burrows & Sawyer 1998; Hannestad & Raffelt 1998; Buras et al. 2003; Itoh et al. 2004; Marek et al. 2005; 
Langanke et al. 2003, 2008 for the extensions to the standard opacities). Neutrino oscillations are presently neglected 
despite the fact that the classical MSW effect and collective neutrino-antineutrino flavor oscillations may possibly have 
an important bearing at least on the observable neutrino properties. However, considering that our understanding of 
the latter phenomenon is by no means complete, and that the classical MSW effect is only relevant far outside the 
supernova core, it is not sure whether neutrino oscillations are dynamically relevant for the supernova problem at all. 
Moreover, no self-consistent dynamical simulations including these effects have been carried out so far; recent studies 
have contented themselves with post-processing progenitor or simulation data (Duan et al. 2008; Lunardini, Miiller, & 
Janka 2008). Thus, our treatment of neutrino physics, while certainly not complete, is state-of-the-art for dynamical 
supernova simulations, and the same is true for the different high-density EoSs (Lattimer & Swesty 1991; Shen et al. 
1998; Hillebrandt & Wolff 1985) available in our code. 

Handling the neutrino transport with the "ray-by-ray-plus" scheme of Buras et al. (2006b) introduces approximations 
on two different levels: First, the variable Eddington factor technique is not strictly equivalent to full Boltzmann 
transport even in spherical symmetry, since only a simplified "model" Boltzmann eqiiation is solved to provide the 
closure for the moment equations. This shortcoming is mitigated by the fact that comparisons with analytic solutions 
(Rampp & Janka 2002) and with Boltzmann solvers based on the 5Ar-method (Liebendorfer et al. 2005) show that the 
approximate solution is entirely adequate in many cases. Second, the "ray-by-ray-plus" approach ignores the lateral 
propagation (but not the lateral advection) of neutrinos. However, the analysis of Buras et al. (2006b) at least suggests 
that this approximation does not affect the heating in the gain layer dramatically, and one may therefore expect that 
it captures the overall dynamics during the proto-neutron star accretion phase and also during the subsequent wind 
phase quite adequately. Unfortunately, no detailed comparison of the "ray-by-ray-plus" method with true multi- 
angle transport or alternative approaches, like multi-group flux-limited diffusion (MGFLD), which would allow us to 
better assess the quality of our method, is available as yet. So far, the different approaches to multi-dimensional 
transport hitherto pursued by supernova modellers all have their individual strengths and weaknesses: The only true 
multi-angle simulations (Livne et al. 2004; Ott et al. 2008) for example, have been carried out at the expense of 
neglecting inelastic scattering reactions, and MCFLD has been shown to smear out lateral variations of the radiation 
field quite considerably (Ott et al. 2008), while the "ray- by-ray-plus" method probably tends to overestimate such 
variations. Still, simulations with ray-by-ray transport show effects that are qualitatively similar to those seen by 
Ott et al. (2008) with multi-angle transport. Global asymmetries of the neutrino flux are preserved better than with 
MGFLD, for example. At any rate, ray-by-ray transport can be classified as a state-of-the-art method in the field, 
notwithstanding the fact that it is certainly limited to scenarios where the anisotropies in the matter distribution (at 
least out to the neutrinosphere) do not become exceedingly large. Compared to other radiative transfer methods used 
in multi-dimensional relativistic simulations up to now like the deleptonization scheme of Liebendorfer (2005), which 
was first used in GR by Ott et al. (2007b), or the grey two-moment formalism (using the Eddington approximation) 
by Farris et al. (2008) - it is a much more ambitious and accurate approach to neutrino transport in core collapse 
supernovae. 

Given the limitations of the "ray-by-ray-plus" method, the use of the CFG approximation to the gravitational field 

equations is quite natural: CFC deteriorates in accuracy for strongly deformed, rapidly rotating matter configurations, 
where the "ray- by-ray-plus" method is not applicable anyway - there is little to be gained by a more accurate treatment 
of GR in that regime, ft is notcnvorthy, however, that the CFC approximation has been found to be accurate for a 
wide range of (sometimes extreme) astrophysical scenarios, such as rapidly rotating neutron stars (Cook, Shapiro, & 
Teukolsky 1996), rigidly rotating disks (Kley & Schafer 1999), and rotational core collapse (Shibata & Sekiguchi 2004; 
Cerda-Duran et al. 2005; Ott et al. 2007a, 2007b; Dimmelmeier et al. 2007). In the core collapse scenario, CFC still 
appears to be adequate (Dimmelmeier et al. 2006) for the extremely rapidly rotating models of Shibata & Sekiguchi 
(2005); even toroidally deformed neutron stars obviously pose little problem to CFC. These facts indicate that in 
typical supernova simulations, where the neutron star rotates moderately at most, CFC will not be the weakest link 
in the chain of approximations, and that, at least for our purpose, GR effects are treated correctly in our code. All in 
all, VERTEX-CoCoNuT can thus be regarded as a truly relativistic radiation hydrodynamics code for core collapse 
supernovae, featuring a sophisticated neutrino transport scheme and an up-to-date treatment of the microphysics. 

3. CODE TESTS IN SPHERICAL SYMMETRY 

3.1. Simple Test Problems - Radiating Spheres 

Stationary solutions for the propagation of radiation from an extended spherical central source into vacuum furnish 
one of the simplest kinds of test problems in radiation hydrodynamics. Although not overly challenging per se, they 
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Fig. 1. — Left: Setup and results for test problem 1: The radial velocity profile used for the test is shown in the upper panel. The middle 
and bottom panels show the luminosity L and mean energy (e) in the moving observer frame as obtained with VERTEX-CoCoNuT 
(crosses), as well as the analytic solution (solid lines). For convenience, L and (e) are given in units of their values Lq and (e)^ just outside 
the radiating sphere. Right: Setup and results for test problem 2: The lapse function a and conformal factor 4> are shown in the upper 
panel. The middle and bottom panels show the redshift-corrected luminosity Lrs and the mean energy (e) in the frame of an Eulerian 
observer as obtained with VERTEX-CoCoNuT (crosses), as well as the analytic solution (solid lines). For convenience, L^b and (e) are 
given in units of their values at infinity Lao and (e)^. 
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Fig. 2. — Redshift-corrected luminosity Lrs and mean energy (e) (both in the moving observer frame) for test problem 3. Crosses indicate 
the results obtained with VERTEX-CoCoNuT, the analytic solutions are drawn as solid lines. For convenience, L^b and (e) are given in 
units of their values at infinity Loo and (e)^^- 

provide a useful consistency clieck for tlie numerical implementation of the relativistic neutrino moment equations. 
The principle underlying these tests may be illustrated by considering the neutrino energy equation ((27| in flat space, 
assuming stationarity and specializing to the laboratory frame (where Vr = 0), 



1 d {r^H^ 

J.2 Qj, 



(55) 



Outside the central source (which may be chosen as a homogeneous isothermal sphere of radius R with vanishing 
scattering opacity and frequency-independent absorption opacity as a very simple case), the collision integral vanishes, 
and hence the lab frame luminosity remains constant, 

icui = 167r^i?eui?'^ = const. (56) 

Likewise, since no energy-exchanging reactions occur, and since all the Doppler terms in the moment equations vanish, 
the mean energy of the radiation is independent of the radius r in the vacuum region. 



(£)c 



HcuifHcui = const. 



(57) 



The luminosity and mean energy as measured by moving observers can be constructed from the transformation 
properties of the radiation moments J, H, and K (see, e.g. Mihalas & Weibel Mihalas, 1984). Since J = H = K ior 
r ^ R, we have „ „ 1 



LW^ 1 



2Vr 
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const, (e) W {I + Vr) — const, 



(58) 



for the luminosity L = IGir^H/c and the mean energy (e) in the frame of a moving observer. 

Test problem 1: We have compared the numerical solution computed with VERTEX-CoCoNuT in the comoving 
frame to this analytical solution for the case of a radiating sphere with i? = 4 km and an observer velocity field with 
a peak value oi v = —0.2c a.t r — 150 km. 



0, 



-0.2c 



r — 135 km 



-0.2c (i 



15 km 
50 km\ 



r < 135 km 

135 km < r < 150 km 

150 km < r 



(59) 



A NEW GENERAL RELATIVISTIC SUPERNOVA CODE 



13 



I I I I Mill — I I I I Mill — I I I 1 1 nil — I I I 1 1 nil — r 



I I I 1 1 Mil — I I I 1 1 nil — I I I 1 1 nil — I I I 1 1 III! — r 



0.45 




P, [g/cm 



P, [g/cm ] 



Fig. 3. — Comparison of the entropy per baryon s (left panel) and the electron and lepton fraction Ye and Yi^p (right panel) versus central 
density pc during collapse in AGILE-BOLTZTRAN (solid lines), VERTEX-CoCoNuT (dashed), and VERTEX-PROMETHEUS (dotted 
and dash-dotted for the effective potentials R and A). A slightly lower central electron and lepton fraction than with AGILE is seen for 
the VERTEX codes. After trapping, the central entropy is conserved to very good accuracy in VERTEX-CoCoNuT. 



This setup mimics the velocity profile typically encountered in the accretion phase of core collapse supernovae with a 
strong shock discontinuity at radii of order 100 km. The left panel of Fig. [1] shows that the numerical solution is in 
excellent agreement with Eq. (j58p . with a maximum deviation of less than 2%. Considering that the finite-difference 
representation of the moment equations has not been specifically tuned to reproduce the analytic result, and that 
the energy grid is rather coarse {Ae/e « 0.28), the accuracy achieved by VERTEX-CoCoNuT seems more than 
adequate, and indicates that the velocity-dependent terms in the moment equations have been correctly implemented. 

Test problem 2: The implementation of the gravitational redshift terms can be tested in a similar fashion: We 
now consider the propagation of radiation from an isothermal sphere in a curved spacetime, assuming that the metric 
is diagonal and isotropic, i.e. ds^ 



-a' 



and number equations (obtained from Eqs. 
directly in the vacuum region. 



(EID and (ESI) 



de-" + sin^ 
after setting 



The frequency-integrated neutrino energy 
= Vr = and W = 1) can then be solved 



b'^nr'^ = const.. 



b'^Hr^ = const. 



a (e) = const. 



(60) 



To test our numerical scheme, we have chosen a lapse function a{r) from a pseudo-relativistic core collapse simulation 
with the VERTEX-PROMETHEUS code (Rampp & Janka 2002) (with a minimum value of a(0) « 0.72), thus 
mimicking a poten tial well of similar strength as in real applications. For convenience, the conformal factor (j) is 
set to (f) = V as in the weak-field limit. Again, the luminosity and mean energy are in very good agreement 
with the analytic results for this test case, as shown in the right panel of Fig. [T] The redshift-corrected luminosity® 
Ly-s — a (l67ra0^ i?r^) , which should be constant outside the radiating sphere for our specific choice for a and cj), is in 
fact constant to within 0.02%. 

Test problem 3: Test problem 1 and 2 can be easily combined to serve as a more stringent check for the implemen- 
tation of the moment equations (as more and more non- vanishing terms enter), if we consider non-Eulerian observers 
in a curved spacetime. The redshift-corrected luminosity and mean energy as measured by an observer moving 
with a velocity iv then obey the following simple equations, 

1 + Vr ^ _ o >aI + Vr 



-Lis — IOttq! 



1 



-Hr — const, Wa (1 -I- Vr) (e) ~ const. 



(61) 



\ — Vr 

With the same choice for Vr-, a, and (/) as in test problem 1 and 2, VERTEX-CoCoNuT once more reproduces the 
analytic solution to very good accuracy (see Fig. [5]). 



3.2. Core Collapse in Spherical Symmetry 
3.2.1. Aims and Test Setup 

While the accuracy of VERTEX-CoCoNuT for the simple test problems presented in the last section is reas- 
suring, only a comparison with existing one-dimensional general relativistic neutrino transport codes can serve as a 
true shakedown test for future applications in "uncharted territory" (i.e. multi-dimensional GR neutrino transport). 
Fortunately, such codes are available (Yamada et al. 1999; Liebendorfer et al. 2001b; Liebendorfer, Mezzacappa, & 
Thielemann 2001a; Liebendorfer et al. 2004). A detailed comparison of the relativistic AGILE-BOLTZTRAN and 
the Newtonian VERTEX-PROMETHEUS code (AGILE and PROMETHEUS for short) has aheady been carried 

* The uncorrected relativistic luminosity in a CFG spacetime, i.e. the energy passing through a sphere around the origin of radius r per 
unit coordinate time is just Idiraij)'^ Hr"^ (note that in the adopted gauge, the circumferential radius is not given by r, but by r(fp). If the 
spacetime is static, gravitational redshift corrections can be taken into account analytically by including another factor a. 
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Fig. 4. — Profiles of the entropy per nucleon s (upper left), the electron fraction Ye (upper right), the radial velocity Vr (bottom left), 
and the density p (bottom right) for model G15 at a time of 3 ms after bounce, obtained with different neutrino transport codes: AGILE- 
BOLTZTRAN (sohd), VERTEX-CoCoNuT (dashed), and VERTEX-PROMETHEUS (case A) (dash-dotted). For the sake of clarity, 
case R is not shown; the profiles are very similar to case A. 
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Fig. 5. — Left panel: Ledoux criterion 20 ms after bounce for the relativistic VERTEX-CoCoNuT run (solid) and for the pseudo- 
relativistic VERTEX-PROMETHEUS (case A) run (dashed). For convenience, Cl is scaled to the local baxyonic mass density p. The 
region between rd = 45 km and r^f = 60 km is convectively unstable in the CoCoNuT run, while no extended unstable region is present in 
the PROMETHEUS run. Right panel: Time evolution of the shock position for model G15 in AGILE-BOLTZTRAN (solid), VERTEX- 

CoCoNuT (dashed), and VERTEX-PROMETHEUS with the effective potentials A and R (dotted and dash-dotted, respectively). 

out by Liebendorfer et al. (2005) and Marek et al. (2006). Since Liebendorfer et al. (2005) provide a useful frame- 
work for assessing the quality of a newly developed code, we have repeated their run G15 with CoCoNuT, and have 
compared our results to those obtained with AGILE. As the neutrino transport module in VERTEX-CoCoNuT is 
largely identical to that in VERTEX-PROMETHEUS (apart from the underlying moment equations), there is an 
additional payoff from this testing strategy: in a number of cases, the origin of the differences between the three codes 
can be pinpointed more accurately with a third code available. In addition to the PROMETHEUS run discussed in 
Liebendorfer et al. (2005), which was carried out using a potential derived from the TOV equation ("case R"), we also 
considered a PROMETHEUS run with a different gravitational potential ("case A") proposed by Marek et al. (2006). 
In total, our analysis therefore encompasses four different simulations^ of model G15 with AGILE, CoCoNuT, and 



® The data from the AGILE and PROMETHEUS (case R) runs have been taken from the online material provided in the electronic 
version of Liebendorfer et al., 2005. 
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Fig. 6.— Neutrino luminosities for model G15 obtained with AGILE-BOLTZTRAN (solid), VERTEX-CoCoNuT (dashed), and 
VERTEX-PROMETHEUS (dotted and dash-dotted for the effective potentials R and A) as measured by an observer at r = 500 km 
in the comoving frame. Electron neutrino luminosities are shown in the upper panels, electron antineutrino and /i/r neutrino luminosities 
in the lower left and lower right panel. 
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Fig. 7. Left: Neutrino rms energies for model G15 obtained with AGILE-BOLTZTRAN (solid), VERTEX-CoCoNuT (dashed), and 
VERTEX-PROMETHEUS (dotted and dash-dotted for the effective potentials R and A), sampled at r = 500 km in the comoving frame. 
The three panels show the rms energies for electron neutrinos, electron antineutrinos and /i/r neutrinos (from left to right). 

PROMETHEUS (with two different choices for the gravitational potential). 

Following Liebendorfcr ct al. (2005) as closely as possible, wc launch our simulation from the 15Mq progenitor 
sl5s7b2 of Wooslcy & Weaver (1995), and use the same set of neutrino interactions rates^" and the same high-density 
EoS of Lattiincr & Swesty (1991). Our treatment of the EoS differs slightly from both AGILE and PROMETHEUS 
for this test run: While using the EoS of Lattimer & Swesty (1991) above a threshold density of 6 x 10^ g cm~^, 
nuclear burning is switched off below 3 x 10^ g cm~^ (as in AGILE), and the composition of the progenitor is 

1" The details of the implementation differ, however: In particular, both VERTEX codes simplify the physics in the /i/r-neutrino sector 
by not treating neutrinos and antineutrinos separately. 
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Fig. 8. — Radial profiles of the velocity v (dashed lines), the density p (solid lines, upper left), the electron fraction Ye (solid lines), the 
entropy s per nucleon (dashed lines, bottom left), the neutrino luminosities L^, (upper right), and the neutrino rms energies (e) (bottom 
right) in AGILE-BOLTZTRAN (black), VERTEX-CoCoNuT (green), and VERTEX-PROMETHEUS with the eff'ective potential A 
(blue) at a post-bounce time of 100 ms. Neutrino luminosities and energies are given in the comoving frame; solid lines are used for electron 
neutrinos, dashed lines for electron antineutrinos, and dash-dotted lines for fi/r neutrinos. 

advected with the fluid (different from AGILE, which considers only one nucleus, ^^Si, in that regime). Between those 
transition densities, the composition is taken from a 17-species NSE table^""^. An Eulerian grid is used both for the 
hydrodynamics and the neutrino transport (400 and 234 zones, respectively), and the energy resolution is the same as 
in PROMETHEUS (17 zones). 

Although the microphysical input is almost identical in all three codes, it is crucial to bear in mind that they differ 
considerably from each other in their approach to neutrino radiation hydrodynamics before proceeding with a com- 
parison of their results: AGILE solves the relativistic Boltzmann equations directly by means of a discrete- angle (Sn) 
method, and relies on an artificial viscosity scheme with an adaptive (moving) grid for the hydrodynamics. VERTEX- 
PROMETHEUS employs a variable Eddington factor technique for the neutrino transport and the Piecewise Parabolic 
Method (PPM) with an exact Riemann solver for evolving the equations of hydrodynamics; it uses a number of approx- 
imations to incorporate general relativistic effects, while remaining essentially a Newtonian code. The key element of 
its GR treatment is a modification of the Newtonian gravitational potential that leads to the same hydrostatic stellar 
structure equation (TOV equation) as in general relativity (VERTEX-PROMETHEUS "case R") or to a sightly 
modified TOV equation that improves the agreement with the relativistic case^^ ("case A"). VERTEX-CoCoNuT, 
while based on the same general approach to neutrino transport and hydrodynamics as VERTEX-PROMETHEUS, 
accounts for general relativistic effects exactly in spherical symmetry (apart from the computation of the Eddington 
factors). However, the approximate Riemann solver used in CoCoNuT (HLLE for this run) is somewhat less accurate 
than the exact solver in PROMETHEUS. Moreover, the two relativistic codes CoCoNuT and AGILE differ in their 
choice of gauge and slicing conditions: Ambiguities arising from this fact must be handled carefully. 

In most instances, the different gauges employed by CoCoNuT and AGILE do no pose any problem because the 
conversion of the gauge-dependent quantities is usually straightforward. Thus, the areal (circumferential) radius rcf , 
which is used by Liebendorfer et al. (2005), can easily be computed from the isotropic radial coordinate r, which 
actually appears as r in the CFC metric: 

rcf = (62) 
We also follow Liebendorfer et al. (2005) in choosing the rate of change of the areal radius per proper time as radial 

This procedure has been adopted to ensure that material in grid cells passing from the high-density to the low-density regime has a 
physically reasonable composition. 

For a heuristic motivation of the modifications in Case A, see Marek et al. (2006). 
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Fig. 9. — Radial profiles of the velocity v (dashed lines), the density p (solid lines, upper left), the electron fraction Ye (solid lines), the 
entropy s per nucleon (dashed lines, bottom left), the neutrino luminosities L^, (upper right), and the neutrino rms energies (e) (bottom 
right) in AGILE-BOLTZTRAN (black), VERTEX-CoCoNuT (green), and VERTEX-PROMETHEUS with the eff'ective potential A 
(blue) at a post-bounce time of 150 ms. Neutrino luminosities and energies are given in the comoving frame; solid lines are used for electron 
neutrinos, dashed lines for electron antineutrinos, and dash-dotted lines for fi/r neutrinos. 

velocity v. However, the identification of time slices for the comparison of radial profiles of density, luminosity, etc. is 
non-trivial. In principle, a complete reconstruction of the entire spacetime would be necessary to allow for a rigorous 
comparison of the code output. Since the data provided in Liebendorfer et al. (2005) are insufficient for such a 
reconstruction, we are forced to identify time slices by the elapsed coordinate time t (i.e. proper time for a non-moving 
observer at infinity). This is actually less critical than it might appear if we synchronize the different codes at bounce 
and use the post-bounce time tph as time coordinate, because the system soon settles down to a quasi-stationary state 
for which the slicing conditions in AGILE and CoCoNuT are identical to very good approximation. 

Fortunately, we can even form an estimate of the asynchronicity that develops between the CoCoNuT and AGILE 
runs, by considering the difference of the elapsed proper time tpropor at the center of the proto-neutron star after a 
given interval St of coordinate time. Since the velocity and shift vector vanish at the center, ^proper is given in terms 
of the central lapse function ac by r 

ipropcr = / acdt. (63) 

The elapsed proper time between the CoCoNuT and AGILE runs therefore differs by 

Aiproper = J (a,C°C°NuT _ ^AGILE) (g4) 

Taking the difference of the central lapse function at bounce (obtained by extrapolating to Vcf = 0) as an upper 
limit, we find that the two runs become asynchronous by at most another 2 ms every 100 ms at the center of the 
proto-neutron star. Since the |q,CoCoNuT _ q,agilE| a,gain becomes smaller after bounce, one may conclude that the 
two runs are never more than a few ms out of sync until the end of our simulation 250 ms after bounce. Backed by 
these findings, we can undertake a detailed comparison of the simulations in CoGoNuT and AGILE (and, of course, 
also in PROMETHEUS) during the collapse, bounce, and post-bounce phases. 

3.2.2. Collapse, Bounce and Neutrino Burst Phase 

Relativistic effects are of minor importance until the late stages of collapse when the density reaches several 
10^'^ g cm^^; only then do the infall velocity v and the compactness parameter M/R of the iron core reach val- 
ues of the order of 0.1 (in relativistic units). Consequently, the three codes are in excellent agreement during most of 
the collapse phase. Slight differences in the central electron and lepton fraction after trapping (i.e. at a density of a 
few 10^^ g cm~'^) can be discerned in Fig. [31 Since the VERTEX-based codes consistently produce somewhat lower 
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values, it is likely that this difference stems from small differences in the treatment of the neutrino interaction rates. 
It should be noted that the new VERTEX-CoCoNuT code conserves the central lepton fraction and entropy to very 
good accuracy, thus passing another important test. 

The good agreement between the different codes persists until shock formation (defined as the instant when the 
specific entropy inside the sonic point first reaches 3fcB/nucleon), which occurs at an enclosed mass of O.SSA/q in 
AGILE, CoCoNuT and PROMETHEUS (case A), and at a slightly smaller mass of O.49M0 in PROMETHEUS 
(case R). Interestingly, important special relativistic effects can be seen shortly after bounce (see Fig.|3]): Within 3 ms 
the shock has propagated to the mass shell m ~ I.OMq. The velocity profile for the best effective potential (case A) of 
Marek et al. (2006) still matches the profiles from AGILE and CoCoNuT perfectly, but the entropy profile left behind 
by the shock differs appreciably. As long as the shock travels through an optically thick medium (to neutrinos), the 
post-shock entropy is higher by up to CSfce in the relativistic case (AGILE, CoCoNuT), indicating that the shock is 
initially stronger. The post-shock velocity and the velocity jump across the shock indeed remain larger in CoCoNuT 
than in PROMETHEUS out until m 0.9Mq - velocity profiles from AGILE during the relevant phase are not 
available from Liebendorfer et al. (2005), unfortunately. Once the shock reaches m « O.75M0, the combination of 
the relativistic jump conditions and the transition to the semi-transparent regime results in the formation of a shallow 
trough in the entropy profile, which is centered around m w O.OA/q. On the other hand, the specific entropy increases 
almost monotonically with m out to the shock radius in that region in the PROMETHEUS run, and the post-shock 
entropy is even a little higher than in the GR runs. 

Although a detailed analysis of the cause of these differences is rather complicated, it is at least possible to account 
for the higher post-shock entropies inside m « 0.75Mq in the relativistic case. In this region, the material is still 
optically thick, and the different jump conditions thus remain the only distinguishing factor from the Newtonian 
case. One can therefore estimate the relativistic correction to the post-shock entropy by comparing the Newtonian 
and relativistic form of the jump conditions and applying a perturbative analysis: The Newtonian Rankine-Hugoniot 
conditions and their relativistic analogue (Taub 1948; Thorne 1973) can be recast into equations for the difference of 
the pre-shock and post-shock pressure (Pi and P2) and specific enthalpy {hi and /12). In the relativistic case, P2 — Pi 
and /i2 — hi are given by, 

hi/pi - /l2/P2 



and 

h2 — hi = 



hi 



hi\^ //i2 ^ ^ 



Pi/ \P2 



(66) 



Here, M is the baryonic mass flux through the discontinuity injdie rest frame of the shock, which is given in terms of 
the three-velocities vi and V2 and the Lorentz factors Wi and W2 on either side by 

M = piWiVi = P1W2V2. (67) 

The three- velocities and Lorentz factors are measured in the rest frame of the shock (which is denoted by the tilde) . 
The Newtonian limit of the jump conditions can then be found by setting /ii — > 1, ft,2 -> 1 in Eqs. (j65|) and (|66|) 
and Wi — ?► 1, W2 — )■ 1 in Eq. (|57)) . By linearizing Eqs. (|65IfH7)) around h = 1 and W = 1 and around a Newtonian 
solution for h and P, one obtains an estimate for the relativistic post-shock enthalpy and pressure (as well as the 
other thermodynamic variables). For the post- and pre-shock conditions in the PROMETHEUS run, this analysis 
predicts post-shock entropies that are higher by up to As = 0.4fcB in the relativistic case, mostly during a very brief 
period where the infall velocity vi is as large as |ui| — 0.34 in the rest frame of shock. This prediction is in rough 
agreement with the actual simulation data. As the shock moves further out beyond m = 0.75Mo, \vi\ decreases and 
the relativistic correction to the post-shock entropy soon becomes negligible (e.g. As < O.OSfce for < 0.2), resulting 
in a negative entropy gradient in the range m sa O.SM© . . . 0.9Mq. For m > 0.9Mo the relativistic post-shock entropy 
then matches the Newtonian one very closely. Although neutrino losses from the post-shock region also start to play 
a role at this stage, this reduction of As with decreasing thus provides a qualitative explanation for the observed 
trough in the entropy profile. 

The effects of relativistic shock propagation are also visible in the density stratification between the proto-neutron 
star and the shock (Fig. 21), which is more shallow in the region outside r^f ~ 25 km (m « 0.7), i.e. just outside the 
region of the first entropy peak in the relativistic models . Such small differences may seem unimportant in the case 
of spherically symmetric problems - in particular when they are washed out after a few tens of ms - but they can be 
of some relevance in multi-dimensional simulations, where they can affect the growth of hydrodynamic instabilities. 
Our ID simulations already give a strong indication that this is indeed the case: Convective instability is expected if 
the Schwarzschild-Ledoux criterion Cl (cp. Schwarzschild, 1958) is positive. Cl is given in terms of the density p, the 
pressure P, the specific internal energy density e, and the speed of sound Cg by 

C,.2££±i)-l|^. (68, 

in the relativistic case (Thorne 1966), or ^ Qp 

dr^f cl dr^i ' ^^^^ 
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in Newtonian hydrodynamics^"^. The left panel of Fig. [5] reveals a convectively unstable layer at r^i ~ 50 km with 
a thickness of more than 10 km at a time of 20 ms after bounce in the CoCoNuT run, which is not present in the 
PROMETHEUS run. 

The stronger deleptonization in the region around m = O.9Af0 also leads to a visibly different evolution of the 
electron neutrino luminosities during the neutronization burst (first panel of Fig. Although the peak luminosity 
ibuist ~ 3.8 X 10^^ erg s~^ in PROMETHEUS (case A and R) agrees well with the one in AGILE, the radiated 
energy between ipb = 1 ms and tpb = 8 ms is smaller by a around 10%. Conversely, the total energy radiated in 
i^e-s in CoCoNuT during the burst is in perfect agreement with AGILE despite the higher peak value of Lburst ~ 
4.3 X 10^^ erg s~^. The different shape of the burst in AGILE and GoCoNuT is probably due to the choice of the 
radial grid, or the different accuracy of the discretization scheme (first order vs. second order in space), which both 
infiuence the numerical diffusivity of the transport scheme and thus lead to a different degree of broadening during the 
propagation of the neutrino burst from the decoupling radius to the observer radius at — 500 km (cp. Liebendorfer 
et al., 2004 for an analysis of the broadening of the burst). The treatment of the model Boltzmann equation in 
CoCoNuT may also be an important factor, either because of the approximations in the tangent ray-scheme (which 
may affect the transition from the optically thick to the optically thin regime) or because of the superior angular 
resolution of the tangent-ray scheme at large radii compared to the ^Ac-method used in AGILE, but these factors 
probably play a minor role. 

3.2.3. Accretion Phase 

The post-bounce evolution of the luminosities of all neutrino flavors (see also Fig. 16]) is again very similar in Co- 
CoNuT and AGILE. During the first 150 ms after bounce the luminosities in CoCoNuT tend to be slightly higher 
than in AGILE. The agreement is about as good as for PROMETHEUS (case A) and significantly better than for 
case R, where the strength of the gravitational potential and hence the accretion luminosity is considerably overes- 
timated. While PROMETHEUS (case A) is always very close to both CoCoNuT and AGILE, it is interesting to 
note that only the two last-named codes show a very abrupt transition from the fast rise in the v^/vt luminosity until 
tpi, = 20 ms to the subsequent plateau phase. The small differences between CoCoNuT and AGILE are well within 
the range expected for different codes (cp. the Newtonian run N13 in Liebendorfer et al., 2005) or even for one code 
at different resolutions (cp. Marek, 2007). As in PROMETHEUS, there is a noticeable drop in the luminosities of 
all flavours (z/g and Pe in particular) between ipb = 150 ms and tph = 200 ms - i.e. at a time, when the shock reaches 
the oxygen-rich silicon shell ~ which is far less pronounced in AGILE due to the different treatment of the nuclear 
composition in the low-density regime and the possible superiority of the Riemann solver methods in PROMETHEUS 
and CoCoNuT in following the composition discontinuity at the Si-SiO interface. 

There is also good agreement between the three codes concerning the spectral properties of the emitted neutrinos 
(see Fig. [T]). The root mean square (rms) energies of electron neutrino and antineutrinos as measured by an observer 
at Tci = 500 km (in the comoving frame) in CoCoNuT differ from those obtained with AGILE by 1 MeV at most. 
For electron neutrinos and antineutrinos the agreement is better than in PROMETHEUS (case R and A), where the 
rms energies are either over- or underestimated. On the other hand, the fi and r neutrinos are a little more energetic 
in CoCoNuT than in AGILE during the first 15 ms after bounce, whereas PROMETHEUS (case A) is very close to 
AGILE. Given the fact that the luminosity of fi and r neutrinos during the first 40 ms after bounce rises much more 
slowly than in AGILE, it is likely that the excellent matching of f^/i^r rms energies between case A and AGILE is 
due to a cancellation of the error introduced by the approximate treatment of general relativity (leading to lower rms 
energies) and differences in the implementation of f^/z^r interaction rates or the transport schemes (leading to higher 
rms energies). Conceivably, the higher peak value in CoCoNuT may originate from the different neutrino transport 
algorithms, finite-difference schemes, and numerical grids used in VERTEX and AGILE, which we already suggested 
as reason for the slightly higher peak luminosity of the neutrino burst. 

Similar luminosities and spectral properties already suggest similar neutrino loss rates in the cooling zone and 
similar heating rates in the gain region, therefore the dynamical evolution of the proto-neutron star and the accretion 
shock in GoCoNuT and AGILE should not be very different either. The right panel of Fig. [5] shows that this is 
indeed the case: Until tpb = 150 ms the accretion shock in CoCoNuT is very close to the one in AGILE; the 
radial deviation corresponds to only one or two grid zones. Interestingly, there is a brief transient period of shock 
stagnation at ipb ~ 7 ms (lasting a few ms) both in AGILE and CoCoNuT, which is absent or at least much less 
pronounced in PROMETHEUS (cases A and R). It is conceivable that this feature is connected to the different 
relativistic and Newtonian form of the jump conditions, and to the slightly higher energy loss in electron neutrinos 
during the burst in the GR simulations. In all simulations, the maximum shock radius of Tcf^max ~ 150 km is 
reached at tph = 60 ms ... 80 ms. Afterwards the shock position in CoCoNuT lies somewhere in between the two 
PROMETHEUS runs untfl tpb ~ 150 ms, and is a little closer to AGILE than either of them. The CoCoNuT and 
PROMETHEUS runs then show a transient phase of shock expansion when the infalling Si-SiO interface reaches the 
accretion front; this feature is absent in AGILE (as explained before) because of the different nuclear composition 
of low-density material. At later times (tpb > 220 ms), the shock in CoCoNuT again recedes to the same radius as 

Note that the Schwarzschild-Ledoux criterion reduces to the simple Schwarzschild criterion Cs = ds/dr for a chemically homogeneous 
fluid. The alternative form Cl = ds/dr^f (dp/ds)p y +dYi^-p/ dr^{ (^p/^Vicp) p (Buras et al. 2006b) can not be applied in the relativistic 
case. However, the relativistic Ledoux criterion can'6% expressed in terms of tfie spatial derivatives of s and Y\^-p in the following manner, 
Cl = ds/dr^i (d(p + pe)/9s)p y^^^ + dYi^p/dr^f {d(p -f pe) / dYi^p) ^ p . 
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Fig. 10.— Left: Time evolution of tlie central density pc for model G8.8 with CoCoNuT (solid line) and PROMETHEUS (dashed line). 
Right: Specific entropy s (solid lines) and electron fraction Ye (dashed lines) in the neutrino-driven wind at a fixed radial coordinate of 
500 km as a function of time for model G8.8 in the CoCoNuT (thick lines) and PROMETHEUS (thin lines) runs. 

the shock in AGILE. The overall shock trajectory in PROMETHEUS (case A, but not for case R) is very similar 
to the one in CoCoNuT, which indicates that the effects of general relativity can still be adequately captured by an 
effective gravitational potential in one-dimensional simulations. 

The very close overall agreement between CoCoNuT, AGILE and PROMETHEUS (case A) and the even better 
agreement between CoCoNuT and AGILE is also evident in Figs. [5] and ^ which show radial profiles of selected 
physical quantities at 100 ms and 150 ms after bounce. There are only minute differences in the velocity profiles 
produced by the three codes; the shock front seems to be resolved a little more sharply in CoCoNuT, but this is 
probably incidental, since the sharpness of the jump changes slightly as the shock slowly moves from zone to zone. The 
effects of genuinely relativistic shock propagation in the early post-bounce phase, which had resulted in a somewhat 
different density stratification in the post-shock region in the GR case, have now been largely washed out. This can 
also be seen in the entropy profiles, which are very similar during this phase, except for the fact that PROMETHEUS 
(case A) seems to give slightly lower entropies in the gain region. Interestingly, we now find the electron fraction behind 
the shock to be somewhat lower in PROMETHEUS (case A), contrary to the situation near shock breakout. The 
luminosity and neutrino rms energy profiles show some visible, but not worrisome differences between the three codes; 
e.g. the Doppler jump at the shock is generally a little larger in CoCoNuT than in AGILE, which is to be expected 
however, since there are small differences in the pre- and post-shock velocities. Considering that the three codes use 
different radial grids (and in the case of AGILE also a different grid in energy space), such small discrepancies seem 
unavoidable and are probably not to be ascribed to coding errors in any of them. 

3.2.4. Long-time Tests 

Unfortunately a direct verification of our new code beyond 250 ms by means of a comparison with AGILE- 
BOLTZTRAN is not easily possible, because data for model G15 are not provided by Liebendorfer et al. (2005) 
beyond that point. However, the mere demonstration that VERTEX-CoCoNuT can stably evolve supernova models 
considerably further would be valuable in its own right, and might also add some insights concerning the reliability 
of the effective potential approach used in PROMETHEUS at very late times. We have therefore also repeated one 
of the longest-running ID simulations carried out with VERTEX-PROMETHEUS, viz., a model ("G8.8") of the 
explosion of the S.SM© progenitor of Nomoto (1984, 1987) with an O-Ne-Mg core. To this end, we have included the 
effect of electron captures on ^'^Ne and ^''Mg (Takahara et al. 1989) and the same approximate burning treatment as 
in Kitaura et al. (2006) in the CoCoNuT code. Different from the G15 run and from published long-time simulations 
of the same progenitor model (Fischer et al. 2009; Hiidepohl et al. 2009), we now use a new nuclear EoS kindly 
provided by J. Lattimer (private communication) in the high-density regime; a more extensive discussion of this new 
EoS in the context of supernova simulations will be given in Hiidepohl et al. (2010, in preparation). 

As for the more massive progenitor used in the G15 run, CoCoNuT and PROMETHEUS give very similar results 
throughout collapse and during the short accretion phase of model G8.8, and we therefore refrain from repeating the 
same detailed analysis of these stages as for model G15. Instead, we avail ourselves of the opportunity to compare a 
full GR treatment and the effective potential approach during the post-explosion phase. Due to the structure of the 
progenitor model, the accretion rate drops rapidly as early as 70 ... 80 ms after bounce when the shock reaches the 
edge of the core, and neutrino heating becomes effective in unbinding the post-shock material, resulting in a weak 
(« 10^° erg) explosion (cp. Kitaura et al., 2006; Janka et al., 2008). Interestingly, an important difference between the 
CoCoNuT and PROA4ETHEUS runs emerges at this point: As the shock propagates through the extremely steep 
density gradient at the edge of the core, it accelerates considerably until it reaches the dilute hydrogen envelope of the 
progenitor whose density in the innermost 10*" km is of the order of only 10~^ . . . 10"^ g cm"'^. In the PROMETHEUS 
run, the maximum post-shock velocity Vps becomes superluminal with a maximum value of Vps = 5.3 c, which clearly 
implies that the propagation of the shock through the outer layers of the core and the hydrogen envelope cannot be 
captured correctly with a Newtonian code. In the CoCoNuT run all velocities correctly remain subluminal with a 
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Fig. 11. — Neutrino luminosities and mean energies for model G8.8 obtained with VERTEX-CoCoNuT (solid lines) and VERTEX- 
PROMETHEUS (case A, dashed lines) as measured at r^f = 500 km in the lab frame. The first three panels show the electron neutrino, 
electron antineutrino and ^/r neutrino luminosities; the right lower panel shows the mean energies for neutrinos of all flavours {ue, Pe, !^;i/t) 
as indicated by arrows. Slightly thicker /thinner lines are used for t'ji/r and Pe, respectively, in order to facilitate their visual distinction. 

maximum value of Vps = 0.9848 c (corresponding to a Lorentz factor of = 5.76). However, the shock becomes 
relativistic only at densities below « 10 g cm~^, and at that point the post-shock temperatures are already far too 
low to allow for nuclear burning. The results of Janka et al. (2008) about the nucleosynthesis conditions behind the 
shock in 0-Ne-Mg supernovae therefore remain valid. 

Apart from the different propagation of the shock in the CoCoNuT and PROMETHEUS runs, the subsequent 
evolution of the newly formed remnant in model G8.8 is remarkably similar in both cases. The contraction of the 
proto-neutron star due to the continuing neutrino losses is almost in perfect agreement as illustrated by the left panel 
of Fig. 1101 which shows the evolution of the central density pc- Even more than 4 s after bounce, pc is only w 3% higher 
in the PROMETHEUS run, and the neutron star radii (not shown) are almost identical as well, with a maximum 
deviation of 2%. Slightly larger differences are observed in the neutrino luminosities (see Fig. Ill|) . which are typically 
a little higher in the pseudo-Newtonian PROMETHEUS run from 0.5 to 2 s after bounce, but the difference to the 
relativistic simulation stays low (« 5%). The mean energies of the emitted neutrinos also agree to within 3% during 
the post-explosion phase, and it is only after more than 2 s that the deviation starts to become larger, reaching a 
maximum of 0.4 MeV for the p/r neutrinos. Considering that the finest energy grid zone has a width of 1.3 MeV, this 
close agreement seems quite remarkable. However, despite the fact that the neutrino luminosities and mean energies, 
as well as the structure of the proto-neutron star remain extremely similar, there is a somewhat bigger discrepancy 
between the two codes concerning the properties of the neutrino-driven wind (see right panel of Fig. I10[) . While the 
electron fraction Yg in the wind differs by less than 0.003 at late times, the entropy s reached asymptotically in the 
wind is roughly 7% lower in the PROMETHEUS run. Since s scales roughly with a small negative power of the 
luminosity (s cx see Qian & Woosley, 1996), it seems unlikely that this is a result of the slightly stronger 

emission of neutrinos in the CoCoNuT run. Furthermore, we find the wind entropy to be lower in PROMETHEUS 
already at early times when the luminosities are still almost in perfect agreement with CoCoNuT. Even if the small 
differences in the mean neutrino energies of the order of less than 3% are also taken into account, the different wind 
entropies still cannot be explained on the basis of the neutrino emission alone: Assuming that the dependence of the 
final entropy on the mean energy of electron antineutrinos (e)p is roughly given by s oc {s)^^^^^ (see again Qian & 
Woosley, 1996), the combined effect of different neutrino energies and luminosities could only account for a relative 
difference of < 2% in s in the most optimistic case. This indicates that the Newtonian treatment of the equations 
of hydrodynamics in PROMETHEUS is probably the cause for the observed discrepancy: Although the effective 
potential used in PROMETHEUS reproduces the structure of static neutron stars very accurately, it is obviously 
somewhat less adequate for stationary (but non-static) flows. However, with a deviation of the specific entropy of the 
wind of less than 10%, the best effective potential (case A) of Marek et al. (2006) can still be considered quite reliable 
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Fig. 12. — Left panel; Violation of ADM energy conservation AEadm as a function of time for model G15. ABadm is shown both in 
absolute values and relative to the initial ADM energy. Right panel: Rate of change of the total neutrino energy inside a given radius R 
(solid) compared to the contributions from source terms and surface fluxes (dashed) at tpb = 230 ms. 

even for the wind phase. 

These resuhs are very reassuring in two different respects: The VERTEX-CoCoNuT code has proved a rehable 
and robust tool even for long-time supernova simulations well into the proto-neutron star cooling phase. In addition, 
the effective potential approach, though unable to reproduce all the details of the relativistic G8.8 run exactly, has 
been shown to be very accurate in a regime where it has hitherto remained untested. 



3.2.5. Energy and Lepton Number Conservation 

Our test runs also provide an opportunity to check the ability of VERTEX-CoCoNuT to conserve total energy 
and lepton number. The usefulness of our improved numerical treatment of the relativistic energy equation presented 
in App. r^and the finite-difference scheme for redshift and Doppler terms in the neutrino moment equation described 
in App. [B]can thus be tested beyond the level of simplified test problems. It should be noted that the very fact that 
the new code allows us to quantify the numerical violation of energy conservation unambiguously is a major advantage 
compared to the pseudo-relativistic approach of VERTEX-PROMETHEUS: For the effective potentials of Marek 
et al. (2006) and Miiller et al. (2008), a conservation for the total energy cannot even be formulated. On the other 
hand, in the framework of the 3 + 1 formalism there is a conservation law for the ADM energy (Landau & Lifschitz 
1997; Misner et al. 1973), which is given by (cp. Eq. (fTU)) ) 

OO TT 27r 

Eadm^ J J j (jy' (^E+^^^^r^smOd^dedr, (70) 



in the case of a CFC spacetime. The accuracy of total energy conservation can therefore be checked by monitoring 
Eadm (taking into account surface fiux terms where necessary) at least in spherical symmetry, where the conformal 
flatness condition does not involve any approximations to the Einstein equations. In the multi-dimensional case, 
the CFC system (fTUl - [T^ is not equivalent to the full Einstein equations, and E'adm is therefore not conserved 
in general. Even though, -Eadm can still be computed and remains a valuable diagnostic quantity for two reasons: 
First, i?ADM should still be conserved to good accuracy as long as the deviations from the conformal flatness remain 
sufficiently small, which is probably the case for the generic supernova problem (cp. Sec. I2.4p . Second, whether the 
non-conservation of -Eadm is indicative of a numerical problem with energy conservation or of some failure of the CFC 
approximation is actually of minor importance; an overly large deviation Ai?ADM from the initial value would be a 
sign of trouble either way. 

However, in the one-dimensional case i?ADM ought to remain constant (neglecting surface fluxes), and Ai^ADM can 
be directly used to measure the numerical energy conservation error. The time evolution of A^^adm is shown in Fig. 1121 
After a very rapid increase by 0.5 x 10^^ erg around bounce, A/?adm gradually grows by another 1.5 x 10^^ erg during 
the next 500 ms after bounce^'*. In two-dimensional simulations to ~ 300 ms after bounce, we also achieve about the 
same accuracy. The significance of these numbers is not immediately clear without further analysis. On the one hand 
the fact that Ai?ADM is of the order of the typical explosion energy may seem worrisome, on the other hand Ai?ADM 
is small compared to the gravitational binding energy. Fortunately, the proper scale of reference for Ai?ADM can be 
worked out by determining the cause of the violation of energy conservation more precisely. As shown in App. |X] for 
the purely hydrodynamic case, |Ai?ADM| should not exceed a few 10^" erg during the bounce phase for our improved 
implementation of the gravitational source term in the energy equation ([7|) , and the additional error accumulated after 
the formation of the proto-neutron star is even smaller. This suggests that the secular increase of Ai?ADM during the 

Interestingly, the evolution of A_Eadm is somewhat different from the relativistic test case in App. [A] where Ai?ADM becomes 
negative at bounce, and afterwards settles down at around —1 X 10^" erg. Imperfect energy conservation in the neutrino sector as described 
in this section may partly account for this. However, the different dynamical evolution of the GR model in App. |A] (where neutrino energy 
losses are neglected completely) may also be of relevance. 
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accretion phase occurs in the neutrino transport sector. To verify this, we speciahze to the case of a static background 
spacetime, where the relativistic evolution equations for the component T"" of the stress-energy tensor T of the matter 
(or the neutrinos) can be written in a fully conservative form by means of the following manipulations, 

(71) 
(72) 

UJb- 

. V ^T-) + A (y^^T") + y^^T- . (73) 




(74) 



d 

Here, T^^" may denote either the stress-energy tensor of the matter or of the neutrinos, and Q is the source term due 
to neutrino-matter interactions, which enters with different signs in both cases (Qmattcr = — Qncutrinos)- The time-like 
components T^^^ of the total stress-energy tensor therefore obey the following conservation law, 

I + A (ay3^r« ) = 0. (75) 

Since our implementation of the gravitational source term in the energy equation ((T]) is essentially equivalent to 
a conservative discretization of Eq. (|74|) . any spurious energy loss or generation in the case of a static background 
spacetime must originate from our finite-difference representation of the neutrino moment equations. This can be 
verified by checking whether the integral version of Eq. (fM| for the neutrinos is fulfilled numerically. More specifically, 
we consider the rate of change of the total neutrino energy inside a given radius i?, which ought to be given by 





gr>^' dr = - [a^r>'] " + / a^Q dr. (76) 



In other words, the^ate of change of the total neutrino energy a^—gT'^^ dr inside R is given by a surface flux 
term — [a^— gT"-'-] ^ and a volume-integrated source term ayryQdr. By plotting the LHS and RHS of Eq. (|76p 
as functions of the radius R as in Fig. 1121 one may determine directly where energy conservation is actually violated: 
Ideally, the numerically evaluated LHS and RHS (solid and dashed lines in Fig. [T^ should be identical. Inside the 
shock at i? « 100 km the agreement is indeed very good, the maximum deviation being 5 x 10^*^ erg s~^. Numerical 
energy conservation is violated more strongly at the shock and in the pre-shock region, and the accumulated error 
amounts to approximately 2 x 10^^ erg s^^, which corresponds nicely to the rate of change of the ADM energy at 
tpb = 230 ms. 

The underlying reason for the numerical violation of neutrino energy conservation can also be worked out easily: 
When deriving Eq. (j74p from the frequency-integrated moment equations in the comoving frame on a static spacetime 
background, all the Doppler and gravitational redshift terms are cancelled by corresponding terms arising from the 
Lorentz transformation to the lab frame and Eq. ((74)) . For the finite-differenced moment equations, these cancellations 
are no longer exact. Further analysis shows that the most critical terms in this respect are the velocity-dependent 
terms of 0{v), which explains why energy conservation is violated at the shock and the pre-shock region, where the 
absolute velocities are highest. In this region, however, neutrino heating and cooling are negligible (except for the 
early post-bounce phase, where there may be some pre- heating), and the dynamically relevant error in the heating 
rate resulting from energy non-conservation in the neutrino sector is therefore much smaller. Thus, the dynamical 
evolution of the accreting proto-neutron star and its environment should not be significantly affected by the secular 
increase of the ADM energy. Furthermore, since Eq. ([76| is violated mainly in the optically thin region, it is quite 
natural to interpret the numerical increase of Eabm of the order of a few 10^^ erg as a luminosity error of a few 
percent and not as spurious energy input. Considering that uncertainties in the microphysics and the discretization 
error probably limit the accuracy of the neutrino luminosities to a similar level anyway, the conservation error observed 
in simulation G15 does not seem very disturbing. Moreover, the problem with the finite-difference representation of 
the Doppler terms in the moment equations becomes less acute in the post-explosion phase once the shock has left the 
grid: In model G8.8, the spurious change of the ADM energy between = 1 s and tpb = 4 s is only 3.7 x 10^^ erg. 
However, since the velocities in the neutrino-driven wind can still reach a few percent of the speed of light, the energy 
error still does not vanish completely. 

Somewhat surprisingly, lepton number conservation is also not fulfilled exactly in models G15 and G8.8; there is a 
secular drift of the order of 1% . . . 2% of the net lepton number luminosity. At first sight, this seems puzzling, since 
the treatment of Doppler and gravitational redshift terms described in App. IB] allows us to solve the neutrino energy 
equation (|27|) in such a way that lepton number is conserved numerically. However, this apparent contradiction can 
be easily resolved: In order to simplify the structure of the Jacobian of the discretized moment equations, we compute 
the flux terms d/de{. . .) using the fiux factor fn from the solution of the model Boltzmann equation and the zeroth 
moment of the radiation intensity J, while the first moment H is used directly in the corresponding terms outside the 
derivative d/de. If fnJ — H held, the cancellation of these terms described in App. [B] would still work, and lepton 
number conservation would be guaranteed; but we find fuJ ^ H for two reasons: First, our Boltzmann equation is not 
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fully consistent with the moment equations (j27p and (j28p in the sense that some terms of the underlying Boltzmann 
equation ((26|) are omitted, and hence the flux factors computed from the solutions of the model Boltzmann and moment 
equations do not agree perfectly, particularly in the vicinity of the shock. Second, fuJ and H are defined on cell 
centers and cell interfaces, respectively, so that the finite interpolation accuracy will, in general, lead to fuJ 7^ H. 
It should be pointed out, however, that such problems are far less severe in the optically thick regime, because there 
H 0. These arguments suggest that the observed imperfect conservation of lepton number conservation is again 
essentially equivalent to a small error in the neutrino (number) luminosity. 

4. CONCLUSIONS 

We have presented a new method for general relativistic neutrino radiation hydrodynamics in core collapse super- 
novae that combines multi-group neutrino transport, an elaborate treatment of the microphysics, multi-dimensional 
GR hydrodynamics and the conformal flatness (CFC) approximation for the gravitational field equations. Our ap- 
proach is a generalization of the variable Eddington factor technique as implemented by Rampp & Janka (2002), 
and relies on the ray-by-ray-plus approach (Buras et al. 2006b) for treating the multi-dimensional case. As in the 
Newtonian approximation (Rampp & Janka 2002; Buras et al. 2006b), the equations of neutrino transport in spherical 
symmetry are solved on different angular bins ( "rays" ) independently, assuming a radial neutrino flux and a parametric 
dependence of the matter and spacetime background on latitude. However, terms for lateral neutrino advection and 
compressional heating of the neutrino gas are also included in the correct relativistic form, as is the acceleration of 
the matter due to lateral neutrino pressure gradients. Although we have considered only the axisymmetric case in this 
paper, the method can be generalized to 3D in a straightforward manner. 

Our approach involves a number of approximations to the full six-dimensional relativistic neutrino transport problem 
on different levels. While the variable Eddington factor technique is well tested in spherical symmetry in the context of 
core collapse supernovae (Rampp & Janka 2002; Liebendorfer et al. 2004), no direct comparisons of the ray-by-ray-plus 
method with multi-angle transport in 2D or 3D are available as yet. However, there is at least qualitative agreement 
with multi-angle transport (see Ott et al., 2008 for results obtained with such an approach and Marek & Janka, 2009 
for a discussion of similarities/dissimilarities to ray-by-ray transport), even in situations where alternative approaches 
like multi-dimensional flux-limited diffusion suffer from excessive lateral smearing of the radiation field. Although 
ray-by-ray transport may not be the method of choice for extremely aspherical configurations (such as accretion tori 
around black holes) , it is probably adequate for the typical core collapse scenario where the proto- neutron star rotates 
only at moderate speeds and is therefore not overly aspherical. As far as the GR field equations are concerned, the 
introduction of CFC is well justified, since it is an excellent approximation even for rapidly spinning compact objects 
for which the ray-by-ray approach would break down. We also note that the relativistic ray-by-ray method could 
easily be generalized to the case where the Einstein equations are solved in spherical polar coordinates without any 
approximations . 

We have implemented our relativistic ray-by-ray-plus method using the framework of two existing codes, the time- 
explicit relativistic hydrodynamics code CoCoNuT (Dimmelmeier et al. 2002a) and the implicit neutrino transport 
code VERTEX (Rampp & Janka 2002; Buras et al. 2006b). Re-utilizing most of the routines of VERTEX allows us to 
include the same up-to-date and well-tested treatment of the neutrino microphysics as in the recent pseudo-Newtonian 
simulations of our group (Buras et al. 2006b, 2006a; Kitaura et al. 2006; Marek, Janka, & Miiller 2009; Marek & 
Janka 2009; Hiidepohl et al. 2009), which greatly facilitates comparisons between the Newtonian and relativistic cases. 

A few improvements were made to the CoCoNuT code, such as the implementation of the approximate Riemann 
solver HLLC (Mignone & Bodo 2005). More significantly, an improved scheme for the energy equation was developed, 
which greatly reduces the numerical violation of total energy conservation for self-gravitating systems. The new scheme 
has also been added to the most recent version of the PROMETHEUS code, and may be of use for other grid-based 
Eulcrian finite-volume codes as well. We also devised a more efficient and consistent treatment for the advection of 
neutrinos in energy space in the VERTEX module. The new method eliminates the need to solve both the neutrino 
energy and number equations to maintain energy and lepton number conservation. 

In this paper, we have limited ourselves to code tests in spherical symmetry before moving on to the first multi- 
dimensional applications of VERTEX-CoCoNuT. While the analytically tractable radiation transport problems 
presented in Sec. 13.11 merely serve as a consistency check for the implementation of the neutrino moment equations, 
we also conducted a full simulation ("G15 run") of the collapse, bounce and post-bounce phases for the progenitor 
model sl5s7b2 of Woosley & Weaver (1995), and compared our results to the relativistic ID Boltzmann solver AGILE- 
BOLTZTRAN (Liebendorfer et al. 2004), building on the detailed comparative study of Liebendorfer et al. (2005). 
Since our code is largely an adapted version of the Newtonian VERTEX-PROMETHEUS code of Rampp & Janka 
(2002), this also allows us to better assess the reliability of the effective potential approach for GR effects (Rampp & 
Janka, 2002; Marek et al., 2006; now also used e.g. by Messer et al., 2007). For this reason, we have also included 
PROMETHEUS runs with two different gravitational potentials (case R and case A) in our comparison. Moreover, 
we performed long-time proto-neutron star cooling simulations of the post-explosion remnant of the S.SAf© progenitor 
model of Nomoto (1984, 1987) with VERTEX-CoCoNuT and PROMETHEUS ("G8.8 run"). 

The results of our tests may be summarized as follows: On the one hand, VERTEX-CoCoNuT is in very good 
agreement with analytic solutions to simple test problems and also with AGILE in the context of a fully-fiedged ID 
simulation from core collapse well into the accretion phase. The agreement with AGILE is slightly better than for the 
PROMETHEUS code with the best gravitational potential currently available. On the other hand, our simulations 
of the 8.8Mq model of Nomoto (1984, 1987) show that PROMETHEUS stih produces results extremely similar to 
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CoCoNuT well beyond the point up to which a direct comparison to the full GR case was hitherto available in the 
work of Liebendorfer et al. (2005). 

The conclusions that can be drawn from these tests are twofold: First, the VERTEX-CoCoNuT emerges as a 
reliable neutrino transport code. Its accuracy has been demonstrated by a thorough comparison of results with a 
well established relativistic ID neutrino transport code. The G8.8 run, covering more than 4 seconds of the proto- 
neutron star cooling phase, also shows it to be sufficiently robust for long-time simulations. For an explicit hydro code 
with an implicit transport solver, the ability to stably evolve a quasi-stationary proto-neutron star over thousands of 
dynamical time-scales and to maintain total energy conservation adequately over such long evolution times cannot be 
taken for granted, but our simulations with CoCoNuT prove that these requirements can be met with an up-to-date 
higher-order Eulerian scheme for GR hydrodynamics. Naturally, the accuracy and robustness of CoCoNuT in the ID 
case will also serve as a point of reference for future multi-dimensional numerical studies of core collapse supernovae 
in general relativity. 

The very good agreement of the rcsuhs of CoCoNuT and PROMETHEUS, both for the G15 and G8.8 run, is 
also noteworthy because it provides hirther evidence that, at least in ID, the evolution of proto-neutron stars can 
be adequately captured by approximating GR effects in a pseudo-Newtonian framework by means of an effective 
gravitational potential. The usefulness and accuracy of the best effective potential of Marek et al. (2006) has now also 
been established for the proto-neutron star cooling phase, and for much longer evolution times and significantly more 
compact neutron stars than in Liebendorfer et al. (2005) and Marek et al. (2006). Despite the smaller differences to 
the GR case seen during the early post-bounce phase, this also gives further credence to 2D simulations using effective 
potentials (Buras et al. 2006b, 2006a; Marek et al. 2009; Marek & Janka 2009; Bruenn et al. 2006; Messer et al. 
2007). 

However, with the advent of VERTEX-CoCoNuT, the choice between rather crude approximations for GR effects 
or a drastic simplification of the microphysics in multi-dimensional supernova studies is now obsolete. The role of GR 
effects in the multi-dimensional case can now be investigated directly with up-to-date microphysics and a sophisticated 
neutrino transport scheme. There are several prominent topics in supernova physics where these new capabilities of 
our code may prove highly valuable, and which we plan to address in follow-up papers. In particular, we intend to 
reinvestigate the dynamics of rotational core collapse and the development of multi-dimensional hydrodynamic insta- 
bilities during the post-bounce phase. The determination of accurate gravitational wave signals - which is connected 
to these issues - is obviously another important goal as well. 

APPENDIX 

A. ENERGY CONSERVATION IN NEWTONIAN AND RELATIVISTIC HYDRODYNAMICS 

A.l. Importance of Energy Conservation in the Supernova Problem 

Conservative discretization schemes for the Euler equations of gas dynamics are one of the cornerstones in modern 
computational astrophysics, both because they have appealing mathematical properties (LeVeque 1998), and because 
they conserve important integral quantities (total rest mass, total energy, etc.) to machine precision. Unfortunately, 
the standard formulation of the Euler equations for self-gravitating fluids is not strictly conservative despite the fact 
that there exists an integral conservation law for the total kinetic, internal, and potential energy for any localized^^ 
matter distribution. 

In the supernova problem, the non-conservative formulation of gravitational source terms in the Euler equations 
raises some critical questions: Can a secular drift of the total energy on the order of 10"''^ erg over a few 100 ms (which 
is a typical value for many simulations; see Liebendorfer et al., 2004; Burrows et al., 2006a; Murphy & Burrows, 2008) 
have a qualitative influence on the propagation of the shock and the development of the explosion, as the explosion 
energy is typically of the same order? Is there an appreciable effect in the first few seconds of the neutron star 
cooling phase? The former question may probably be answered in the negative, becaiise the accumiilated error in 
the gravitational energy liberated diiring the collapse of the iron core and the contraction of the proto-neutron star 
makes itself felt primarily as a small excess or deficit (on the percent level or below) in the internal energy of the 
compact central object (which contains most of the mass on the grid), and not as a large relative error in the region 
around the shock. On the other hand, a secular drift of the total energy can qualitatively affect neutron star cooling 
in the post-explosion phase: Once spurious generation of energy due to numerical inaccuracies outweighs the slowly 
decreasing neutrino losses, the proto-neutron star does not cool any longer (such an effect was actually observed in 
one of our pseudo- Newtonian simulations). 

In order to alleviate the problems associated with the non-conservative standard formulation of gravitational source 
terms in the energy equation, we have developed an alternative scheme for treating the energy equation that significantly 
decreases the secular drift of the total energy. The new scheme can easily be integrated into any Eulerian finite- volume 
method, and is applicable (with only small variations) both in Newtonian and relativistic hydrodynamics. In the 
Newtonian case, the new scheme also has clear advantages when working with effective gravitational potentials that 
do not obey the Poisson equation, despite the fact that a conservation law for the total energy cannot be formulated 
in this case: Using a modified gravitational potential, it turns out that the energy released during the formation of the 
neutron star is in good agreement with the binding energy computed from the Tolman-Oppenheimer-Volkoff equation. 



1^ This restriction is crucial, and has tangible implications: It is not possible, for example, to formulate a conservation law for the energy 
in a general cosmological context (Straumann 2004). 
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at least for neutron stars of moderate compactness (cp. Hiidepohl et al., 2009), while systematic energy conservation 
errors originating from the discretization of the gravitational source term may qualitatively change the outcome of 
numerical simulations. 

A. 2. Energy Equation in Newtonian Hydrodynamics 
A. 2.1. Description of the Algorithm 
Usually the energy equation for a self-gravitating fluid is written as 



de 
dt 



d{e + P) v' 



(9$ 



(Al) 



dx^ dx^ 

Here, e = pv'^ /2 + pe is the total energy density, P is the pressure, is the fluid velocity, p is the (baryonic) mass 
density, e is the specific internal energy, and $ is the Newtonian gravitational potential. Although Eq. (jAll) is not 
written in the form of a pure conservation equation, there is a conservation law for the volume integral of the total 
(i.e. kinetic, internal, and potential) energy (see, e.g. Shu, 1992): In order to obtain this integral conservation law, we 
absorb the source term —pv^d^/dx'^ into the flux term on the LHS, 



and use the continuity equation dp/dt 
term. 
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to eliminate the divergence of the mass flux in the new source 
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AnGp, and the fact that the surface integral, 
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vanishes if the domain of integration is extended to infinity. Eq. (jA4p implies that the total energy is conserved, 
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However, the discretized representation of the total energy, 
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will not be automatically conserved, if the numerical implementation of the gravitational source term is based on 
Eq. (IAl|) . because the discretized source term cannot be manipulated in the same fashion as in the steps from Eq. ()Al|) 
to (jASI) . One way to eliminate any secular drift of the total energy resulting from this, and to avoid the unwanted 
consequences described in Sec. lA.ll would be to use a strictly conservative form of the energy equation (Noonan 1984), 
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Unfortunately, thi^ form of the energy equation has several drawbacks: The "gravitational self-energy" contribution 
l/(87rG') {d^/dx") to the total energy density can exceed the kinetic, internal, and binding energy of the matter by a 
large factor, resulting in considerable round-off errors in these quantities. Furthermore, the flux term contains a mixed 
spatial and temporal derivative, which may introduce numerical inaccuracies. It is also impossible to formulate an 
analogue of Eq. for the effective relativistic potentials commonly used in VERTEX-PROMETHEUS (Marek 
et al. 2006; Miiller et al. 2008) because these do not obey the Poisson equation. 

However, there is a simple and efficient solution to the problem of numerical energy conservatio n th at avoids these 
shortcomings and requires only a minimal modification of the "standard" scheme based on Eq. (IAl|) . The starting 
point is the reformulation of the energy equation in Eq. (jA2l) : Using the continuity equation, the divergence of pv^ can 
be replaced by the time derivative of the density, which is in turn partially absorbed in the time derivative of p^, 



d{e + p<^) d{e + P + p^)v^ 



dt dx^ dt 

The remaining source term now contains the time derivative of the gravitational potential, 
solved with an operator-split scheme in three steps. 



Eq. 



(A9) 

H) can be conveniently 
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Step 1: First, the equations of hydrodynamics are solved with the additional potential energy terms p$, where $ 
itself is kept constant, 

Since /o$ may be considerably larger than the internal energy density e, using e + /5$ as conserved quantity instead of 
e may lead to large round-off errors in e and produce unwanted numerical noise (in particular spurious entropy loss or 
production). To circumvent this problem, the time derivative dp^/dt can be re-expressed in terms of the mass flux 
(as $ = const.), i.e. dp^/dt — ^dp/dt = — <i>V • (pv). 



de d{e + P + p^) dpv^ 



$ = 0, (f> = const. 



(All) 



dt dx^ dx^ 

The value e^^+^Z^) of the energy density after the first fractional step can be obtained by integrating Eq. (jAlip over 
the cell volume V . It is instructive to compare the final formula for the update of e with "standard" implementations 
of the gravitational source term. Working with forward differences in time^^, we now have 



g(«+l/3) ^ g(n) 



[ (e + F) V dA + V / ($e - $fc) pv • dA, 

JdV ,. JAk 



(A12) 



(where Ak denotes the fc-th cell interface, while and $fc are the cell-center and -interface values of $) instead of 



(ti+1/3) 



At 



[ (e + P) V • dA / V$ • pvdV 

Jav J 



(A13) 



for the standard formula tion. Since / {^k — ^c) dA^/V is the component of V$ p ointing into the direction of dAk, 
the update prescription (jA12p can be viewed as an alternative version of Eq. (|A13p , in which the cell-averaged mass 
flux pv has been replaced by an average of the mass flux at cell interfaces. Thus, gravitational energy is only released 
or stored when mass is actually exchanged between different cells. This is the key ingredient for eliminating a secular 
drift of the total energy. 

When the gravitational field is stationary, no further fractional steps are required, and the value of e after one full 
time-step is just e'-"+^-' = e^"^^^^\ If the gravitational potential varies with time, two additional steps are necessary 
once $ has been updated from its old v alue < 1>^"^ to at the end of the time-step. 

Step 2: In order to ensure that Eq. (jAlOp is fulfilled after the update of the gravitational potential $, the energy 
density e^"^^/^' after the second fractional step must be calculated according to 



= (n-H2/3) 



,("-Hl/3) 



$(") _ $("+!) 



Step 3: Finally, the source term pd^/dt has to be accounted for in a third step, 

die + p$) 9$ 



dt 



(A14) 



(A15) 



If centred differences in time are used to obtain second-order accuracy, the energy density at the end of the time-step 

e("+l) ^ ^ 1 (^p(n) ^ p{n+l)^ (^$("+1) _ $(")^ . (Al6) 

Steps 2 and 3 can easily be combined into a compact formula for the update of e after each recalculation of the 
gravitational potential: i , , , , 

g(n+l) ^ g(n+l/3) ^ ^ ^^(n) _ (^$("+1) _ $(") j (Al7) 

A. 2. 2. Condition for the Exact Conservation of Energy 

Under certain conditions, the scheme described in Sec . I A . 2 .Tl conserves the total energy exactly. To demonstrate this 
for the special case of a spherical polar grid, we consider the total energy balance after each sub-step of the algorithm: 
Steps 1 and 2 conserve the (discretized) total volume integral of e + p$, provided that the flux (e + P + p<I>)v vanishes 
at the boundaries of the computational domain. 



(n+2/3) 



_ (n-l-l)^(«+l) 



E 



(An) I (") (()(") \ A 



If e is then updated according to Eq. (|A16|) . we obtain, by combination with Eq. (|A18|) . 



(A18) 



(A19) 



i,j,k 

E 

i,j,k 



(«) 



+ 2 \P{i,j-kr{^-J,k) P(ij,k)^(t.j,k) + Pii,j,k)'^{i,j.k) + P{ij,k)^{i,j.k) ) 



^® First-order forward differences are only adopted for illustration at this point. The implementation in higher-order time integration 
schemes is straightforward. 
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Hence, the total energy after the {n + l)-th time-step is according to the definition in Eq. (jA7l) . 

The sums on the RHS quantify the violation of total energy conservation during one time-step. If they were replaced 
by integrals, we could manipulate them using Green's first identity and the Poisson equation V • V$ — Ai^Gp to obtain 
an exact cancellation: 

j dV^^y" p("+i)$(")dV^ = ^ j (v$(")) • (v$("+i)) '^^"i^ / (^*^"^^^) • (v$^"^) dF = 0. 

Using a special finite-difference representation of the Poisson equation, we can manipulate the sums in Eq. (jA20l) in 
very much the same way as by applying Green's first identity. Let us assume that the source density AirGp is expressed 
as the divergence of the gravitational acceleration g in the following natural manner, 

47rG'pAV(jj^fc) ^g(i+i/2,]M)^A{i+i/2,j,k) - 3(i-i/2j\fc)AA(j_i/2j», 

+ 5(jj + l/2,fe)^^(ij + l/2,fc) ~ 5(jj-l/2,fc)^^(ij-l/2,fe)7 

+ 5(ij,fe+l/2)^^(i,j,fc+l/2) - ff(ij,fe-l/2)^^(ij-,fe-l/2)- (A22) 

Here, A^(i_|_i/2.j,fe) denotes the area of the interface between the cells (i,j,k) and (z -I- \,j,k), and g{i+i/2.j,k) denotes 
the perpendicular component of g thereon; AA(i j_|_i/2,fe) , etc. are defined analogously. Let us further suppose that 
the values of g on any cell interface are obtained by multiplying the difference of $ in the adjacent cells by a time- 
independent factor Qi.j.k, as in the following natural finite-difference representation, 

5(i+i/2j,fc) =Ci+i/2j-,fe (*(i,i,fc) - " ' (^23) 

'2+1 ' i 

5(»j+i/2,fc) =Cij+i/2,fe (*(ij,fc) - =: [Oj+i - 9j) ~ ' (^24) 

5(i,/c+i/2) =Cij,fe+i/2 (*(ij,fe) - *(»,j;fc+i)) sin6>j- (ipk+i - (pk) ^'^^''^'''^ ~ '^{^,j.k+l)) ■ (A25) 

At the outer boundary, we may set g(j\/+i/2 = Cm+i/2 j,fc*(mj,fc) ^{m,j,k)rm/rl^+i/2 to obtain the correct 

monopole moment of the gravitational field. 

The first sum term in the energy balance equation (|A20p can now be written as 

E(n) A -rr 1 \^ / (n) . . (n) « . \ 

P{z,j,k)'^i^U,k)^^h3,k = 2^ [9(^+l/2.J,k)^'^{^+l/2,j.k) ' 9l^-l/2.J,k)^^(^~l/2.J,k) ) <P(ij,fe) 

+ E (4j + l/2,fc)^A^J + l/2.fe) -5(3-l/2,fc)^%J-l/2,fc)) 

+ E (.9(2.^+i/2)^^(^.^fc+V2) - 5£ife-i/2)A^(.,,.-i/2)) (A26) 

The sums on the RHS of Eq. (jA26l) can be rearranged via summation by parts using the finite-difference representation 
for g in Eqs. (|A23HA25l) . For the sum containing the radial components gi+i/2,j,k of the gravitational acceleration g 
on cell interfaces, we have (bearing in mind that A(i/2j,fe) = 0, as i = 1/2 corresponds to the origin): 

E( n^'^'i A A A 4 \ d>("+i) 

(^5'(j+l/2j\fc)^^(»-Hl/2J,fe) - 5(j_i/2j;fc) (*-l/2j:*:) ) *(j,i,fc) 

i:'j,k 

M NO M N O 

= E E E 3(»+l/2j,fc)^^(»+l/2j,fc)*("tfc) ^ E E E 4"-l/2j\fc)^A«-l/2j.fc)*("lfe) 
1=1 j=l fc=l 1=1 j=l fe=l 

M N O M-1 N O 

=EEE5Si/2,..)A^(^+i/2.,.)'j>S:S - E EE5!:^i/2.,.)A^(^+i/2,..)'j>!r^i.) (m 

2=1 j=l fc=l i—1 j — 1 k—1 

M-1 NO NO 

E\~^ (") A A /^r^C'l+l) ^("+1) A I (") A A 

Z^Z^%+l/2,j,fc)^^(*+l/2J,fe) [^(i.j,k) ~ *(i+l,j,/c) ) + Z^Z^%f+l/2J,fe)^^(M-H/2j,fc)1'(Afj,fc) 

i=l j = l fc=l j=l k=l 

M-1 N O ^ , N O 

A A. , . 

/2J,fc)- 



EY^Y^ (") A /I ("+!) 1 Y^Y^ (") A A ("+!) 

2^ 2^ 7 %+l/2,j,fe)^^(«+l/2j,fe)ff(j+l/2,j.fc) + 7 2^ Z^5(M+l/2,j,fc)^^(A^+l/2,j,fc)5'(M+l 

.i=l j = l fc=l ''«+l/2j:fc <,A/+l/2j,A; ^.^^ ^^^^-^ 

The remaining terms on the RHS of Eq. (jA26|) can be manipulated in a similar fashion^'', and we finally obtain an 

-"^^ Note that in spherical polar coordinates further "surface contributions" vanish since j/2,fc) = A(i,JV+i/2,fc) = Oi ^(i,j,0) = ^(i,j,0)> 
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expression that is completely symmetric with respect to the time indices n and n + 1: 



El 



ill ^ o 



(n+1) 

/2j-,fe) 



/2,k) 



^ \ " \ " \ " (n) A A ("+!) 

'7~~~r7Z^ 2^ z^%j+i/2,fe)^^(^j+i/2.fc)%,j+i 

',Jj + l/2,fe ^.^^ 

M N O-l 

'7~7~Z H f ("J^fc+l/2)^^(^J,fc+l/2)5(^,J^fe+l/2) 



w o 



' \ " \ " (") A A 

+ t,f^i/o ^ ^ 5(Af+i/2J,fe)^^(m+l/2j,fc)3(™+i/2,j,fc) 
',Af+l/2j,fe .^-^ ^.^^ 



(A28) 



In effect, we have obtained a finite-difference analogue of Green's first identity with non-vanishing surface contributions. 
The second sum in Eq. (|A20p can be rearranged in precisely the same way and thus cancels the first sum exactly. 
Hence, the total internal, kinetic and gravitational energy is conserved, 

= Ei:l (A29) 



(n+l) 
tot 



A. 3. Energy Equation in General Relativity 

In GR, the situation is somewhat different from the Newtonian case: At first sight, the fact that no sources appear in 
the local conservation law ViyT'"^ = for the stress-energy tensor T^'^ might suggest that the energy and momentum 
equations can be formulated in a strictly flux-conservative form (i.e. without sources). However, an integral conservation 
law can only be formulated (via Gauss' theorem) for divergence-free vector fields on a differentiable and orientable 
manifold (Wald 1984; Straumann 2004); such fields can in general be constructed from the energy momentum tensor 
only if there exists a Killing vector field^*. On the other hand, a conservation law can be formulated for the sum of 
T'^'' and the Landau-Lifshitz pseudo-tensor t'^'^ of the gravitational field (Landau & Lifschitz 1997; Straumann 2004), 



d{-g) {Tf" +tf') 



0. 



(A30) 



Unfortunately, the practical use of Eq. (jA30p is limited because i is a complicated (and non-unique) function of 
the metric and its derivatives. Furthermore, it is impossible to interpret P'''^ as the local energy-momentum vector 
density of the gravitational field due to its non-uniqueness and non-tensorial character; only the integral -Eadm — 
J (r^'^ -I- t^") dV^ (known as ADM energy) over t he en tire spatial domain in an asymptotically flat spacetime has 
a well-defined physical meaning. In general, Eq. (jA30P will only be useful for special gauge choices and symmetry 
assumptions, where the resulting source- free energy equation is not overly complicated (see e.g. Romero et al., 1996; 
Liebendorfer et al., 2001a). The total volume integral Eabm (which can also be expressed in different ways) is, of 
course, an important diagnostic quantity for assessing the quality of energy conservation in a numerical code. 

Because of these obstacles we construct an improved scheme for the energy equation in a similar way as in the 
Newtonian case and do not attempt to eliminate source terms from the energy equation altogether. We start from the 
formulation used by Dimmelmeier et al. (2008), 



T 



00 



^ dx ^ ) 



(A31) 



For convenience, we add the continuity equation, and subsume most of the source terms on the RHS under the variable 
Q to arrive at an equation very similar to Eq. (jAll) . 



Di}^ + Pv') 



da 



(A32) 



dt dx"^ 

The first source term on the RHS corresponds closely to the Newtonian source term pv^ / dx"^ . It is the product of the 
total energy flux (rest mass flux in Newtonian gravity) and the flat-space gradient of t he lap se function (corresponding 
to the Newtonian potential). This source term can be eliminated by multiplying Eq. (jA32| by a, 



and $ 



{i,j,0 + l) 




D) 



da 

dx^ 



(A33) 



Special cases comprise: i) the flat Minkowski spacetime with ten Killing fields (implying the conservation of energy, momentum and 
angular momentum), ii) stationary spacetimes with a time-like Killing field (implying energy conservation), and iii) axisymmetric spacetimes 
with a space-like Killing field that creates an SO{2) isometry group (implying the conservation of one angular momentum component). 
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and by pushing the lapse function into the space and time derivatives, 
d^a (r + D) dy^a {tv' + Dv^ + Pv') 



da 



aQ. 



(A34) 



dt dx' ^ ' dt 

Again, the new source term (t + D) da/dt corresponds closely to pd^/dt in the Newtonian case. Since it is numerically 
advantageous to separate the baryonic mass contributions to the total energy (due to the reduction of round-off errors), 
we subtract the continuity equation, and finally arrive at the following alternative energy equation: 

(j (J (ja 

— [^a {t + D)- ^D] + — [V^ {are + aDff ~ Dv' + aPv')] = {t + D) — + aQ. (A35) 



LV / V ■ / V / J ' LV .-^ V ;j V ; V ■ / 

The numerical solution within a finite- volume scheme can again be split into three sub-steps, the first of which 
consists in solving the energy equation (including only aQ as source term) while keeping the lapse function a fixed, 
d d 

— [y/jar + ^{a-l)D] + — [y/^ {arv' + aDv' - Dv' + aPv')] = aQ. (A36) 



dt dx^ 
In order to avoid working with ^/^ar + ^ (a — 1) as conserved quantity, we note that the change in ^ {a — 1) D 
inside a cell volume V during the time-step At can be expressed in terms of the integral of the flux Dv'^ over the cell 
boundary A = dV. The finite- volume version of the energy equation for this cell can therefore be written as 



Pv') dA 



[ aV7T^"+^/^MF= ( a^^T^") dP^ At - ( ^a^ {t 

Jv Jv V JdV 

- I ^a{a-l)DffdA+{a-l) [ ay^Dv'dA+ j aQdV 

JdV JdV Jv 



(A37) 



The surface integrals in the second line of Eq. (jA37p can be combined into a single term containing the difference 
between the lapse function Uc at the cell center and the lapse function ak at the k-th cell interface, 

(A38) 



f ^ f ^r'-"'^dv + At [ QdV 

Jv JV Jv 



'V 

At 
a 



/ ^a^ {tv' + Pv') + V / v^ttfe (ac - ak) dA 
JdV ,. Ja^ 



This formulation requires only a minimal modification of the scheme for the unmodified energy equation (jA31|) . 
In a second step, the lapse function is updated, while the conserved quantity [ar + (a — I) D] is kept constant, 



d_ 

dt 



aT+{a-l)D 



= 0. 



(A39) 



Note that it is convenient to work with the densitized quantities f = ^J^t and D = here because may be 

updated at the same time as a. The discretized version of Eq. ()A39|) . 



has the solution 



^(n+l)^(n+2/3) ^ |^^(„+1) _ ^ ^(„)^(n+l/3) ^ ^^(n) _ ^(„+l)^ 

^ ^^(n+1/3) _^ (^^(n+1/3) ^ _ 



(A40) 
(A41) 



In a third operator-split step, the source term containing ^/J {t + D) da/dt in Eq. (jA35P is taken care of, 

|[V7«r-f V7(«-l)^] =V7(r + i?)^. (A42) 

As in the Newtonian case, the second and third step can be merged, and if centred differences in time are used, r has 
to be updated according to 



7-(«+l/3) 



r(n+l/3) ^ f (n) £)("+!) + £)(") oi^n) _ ^(n-l) 
2 7^) ' 



(A43) 



A. 4. Numerical Tests 

To illustrate the virtues of the new scheme for the energy equation, we have conducted several numerical tests 
in spherical symmetry with different codes. Newtonian core collapse simulations have been carried out with the 
Newtonian version of CoCoNuT and with a self-written second-order accurate hydrodynamics code employing TVD 
reconstruction and the Kurganov-Tadmor central scheme (Kurganov & Tadmor 2000). 

For the CoCoNuT run, we employ the deleptonization scheme of Liebendorfer (2005), and take neutrino energy 
losses into account in the total energy budget, i.e. we subtract the (negative) time- and volume- integrated neutrino 
energy source term from the total energy measured at any given time. The left upper panel of Fig. [13] shows that the 
new implementation still produces a residual violation of energy conservation of about 2 x 10^^ erg around bounce. 
However, the total energy changes by less than 1.7 x 10^^ erg after tpb = 10 ms, which is a considerable improvement 
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Fig. 13. — Total energy _Etot, central density pc (panel c), and total kinetic energy i?kin,ps inside the shock radius (panels b and d, on a 
linear and logarithmic scale, respectively) for a core collapse runs with CoCoNuT in Newtonian gravity, which has been carried out using 
the deleptonization scheme of Liebendorfer (2005), and the progenitor model s20.0 of Woosley, Heger, & Weaver (2002); the energy lost 
in the form of neutrinos is added to the total energy budget. Results obtained with the "standard" and improved scheme for the energy 
equation are drawn as solid and dashed lines, respectively. 
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Fig. 14. — Left panel; Violation of total energy conservation for the collapse simulation of a 4/3-polytrope with the TVD test code. 
Initially the total energy of the polytropc is —1.08 X 10^^ erg. Right panel: Violation of ADM energy conservation for a general relativistic 
core collapse run with CoCoNuT with the "standard" scheme (solid) and the improved scheme (dashed) for the relativistic energy equation. 
The run has been carried out using the parameterized Ye{p) trajectory of Liebendorfer (2005), and the progenitor model s20.0 of Woosley 
et al. (2002). Neutrino energy losses are neglected during the run. 

compared to the old implementation, where the accumulated error exceeds the absolute value of the initial binding 
energy (w 5.8 x 10^*^ erg) within a few tens of ms after bounce. Nonetheless, the reduction of the numerical error is of 
little relevance dynamically: The proto-neutron settles down at a slightly higher central density, which is consistent 
with the smaller total energy obtained with the new scheme. The total kinetic energy in the post-shock region (right 
panels in Fig. [T^ is essentially identical in both cases as well. We also note that by evolving e instead of e + 
round-off errors in the internal energy and the specific entropy can indeed be kept small; we did not find any unphysical 
entropy production during collapse. 

Although a change of the total energy by around 0.3% as in the CoCoNuT simulation is probably more than 
satisfactory for core collapse simulations, this level of accuracy does not exha ust the full potential of our new scheme, 
which can in principle conserve the total energy exactly, as explained in Sec. IA.2.21 In order to verify this, we have 
simulated the collapse of a 4/3-polytrope with our self- written TVD test code using the hybrid EoS of Janka et al. 
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Fig. 15. — Spectral energy flux of electron neutrinos (solid), electron antincutrinos (dashed) and /j/r (dash-dotted) neutrinos from the 
simulation of a cooling neutron star with VERTEX-PROMETHEUS. All quantities are measured at r^f = 400 km in the comoving frame. 
The left panel shows spectra from 1 s after bounce, obtained with the Doppler treatment of Rampp & Janka (2002). The spectra in the 
right panel (from tpj, = 1.5 s) were obtained with our new conservative Doppler scheme using "harmonic" (a = 1/2) interpolation. Aside 
from some secular changes of the electron neutrino spectrum, the main difference is the elimination of the dip in the Ue- and Wfi/ur spectra 
at the neutrino energy £ 80 MeV. 

( 1993) . Different from CoCoNuT, the gravitational potential and acceleration are calculated according to Eqs. (jA22l 
- IA25|) . hence the condition for exact energy conservation is satisfied. Fig. [14] (left panel) shows that the total energy 
changes by less than 10^^ erg (!), which corresponds to a relative error of less than 10~^^ in terms of the initial energy 
(£;tot,i = -1.08 X 10^1 erg) and less than lO^^^^ in terms of the internal energy at the end of the simulation. Thus 
our improved scheme can indeed preserve the total energy almost to machine precision provided that the gravitational 
potential is determined appropriately. 

In addition to these tests in Newtonian gravity, we have also simulated the collapse of model s20.0 in GR using 
CoCoNuT. We employ the Liebendorfer scheme for deleptonization with one minor modification: While reducing Ye 
according to a ye(p)-pai'anieterization, we neglect energy losses due to neutrino emission, because this simplifies the 
total energy budget considerably^^. Ideally, we would expect the ADM energy E'adm, which can be written as 



in a CFG spacetime, to be conserved in such a situation. Fig. [T3] (right panel) shows that neither the "standard" 
scheme nor our new scheme manage to conserve i?ADM perfectly. In both cases, conservation is violated by |Ai?ADM| ~ 
6 X 10^'^ erg around bounce, which is not surprising considering the complicated non-linear nature of the field equation 
and the appearance of gravitational self-energy terms in the sources and in the ADM energy. However, |A£^adm| 
settles down to a value of less than 10^*^ erg within 1 ms after bounce in the simulation using our new scheme, and the 
secular drift is only 10^^ erg for the next 100 ms afterwards. The old scheme, on the other hand, leads to a violation 
of around 9 x 10^° erg and a small but clearly visible drift of Eabm during the post-bounce phase. 

The relativistic test thus confirms our findings in Newtonian gravity: Our reformulation considerably improves 
numerical energy conservation in simulations of self-gravitating systems. The secular drift during the quasi-stationary 
proto-neutron star evolution is not completely eliminated in all cases, but can be easily reduced by an order of 
magnitude, which demonstrates the usefulness of our new scheme particularly for long-time simulations covering many 
dynamical time-scales. 

Setting the neutrino energy source term to zero while retaining a source term for Ye is of course unphysical, but, as far is our test 
problem is concerned, the prescription adopted here is only employed as a means of initiating the collapse. A non-zero neutrino energy 
source term would complicate the energy budget considerably: A ny up date of the internal energy of the fluid (or, equivalently, the entropy) 
would not only change the terms h and P in the integrand in Eq. IA44I locally, but would also affect the conformal factor <f> and the extrinsic 
curvature Kij non-locally. On the other hand, this problem does not arise with full neutrino transport, because the neutrino source terms 
do not change the total (matter and neutrino) energy-momentum tensor T^'^^ in that case, and therefore do not affect the gravitational 
field. 




(A44) 
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B. SIMULTANEOUS CONSERVATION OF ENERGY AND LEPTON NUMBER 

While the monochromatic moment equations ([27l ISS)) and pO) [3T|) for neutrino energy and number transport are 
analytically equivalent, a flux-conservative discretization of either pair of equations does not necessarily guarantee 
neutrino number (or, in the presence of neutrino-matter interactions, lepton number) and energy conservation simul- 
taneously. The reason for this problem, which is inherent to spectral neutrino transport in the comoving frame, can be 
understood by comparing the analytic and discretized forms of the Doppler and gravitational redshift terms (denoted 
by a subscript "D") in the moment equations. In the neutrino energy equation these terms can be written as 

dJ\ dewJ 

where the advection "velocity" in energy space w is a function of the velocities, metric functions, and the Eddington 
factors Jh = H/ J, fx — K/ J , and /l = L/ J. To obtain the corresponding term in the equation for the monochromatic 
neutrino number density J = J/e, we divide by e and apply the product rule, 

dj\ dwJ wJ wJ dm J dew J 



dt J ^ de e e de de 

The Doppler term in the neutrino number equation is thus a pure flux derivative term and cancels out if the equation 
is integrated over e to obtain the evolution equation for the total neutrino number density. On the other hand, if the 
energy equation is discretized with backward differences in time and a flux-conservative form for the energy advection 
term. 



e»+i/2Wi+i/2^^V/2 - e-*-i/2U'«-i/2-^I!_Y/2 ~ WiJf+^Ae,, (B2) 



At 

the discretized Doppler term for J reads 



-xt / - ^^^+1/2-^+12 - ^--1/2^^1 2 - -^-^r V 

n I n I n 1 

= Wj+i/2ei+i/2Ji" 1/2 ~ '^i-i/2ej-i/2 J^i" 1/2 ~ ^iJ'i — — 
£j+i/2 - e-i „+i £i - ^^i-i/2 „+i 

+ W*+l/2^,+i/2 + ; W,_i/2J,_i/2. 

' (B3) 

Here, Ei and £^+1/2 denote the values of the energy at cell centers and interfaces, respectively, and Ae^ ~ £1+1/2 ~£i-i/2 
is the width of the z-th zone in energy space. The source term in the neutrino number equation thus does not vanish 
in the discretized form of the equations, unless the last three terms in Eq. (|B3P cancel. Different methods have been 
proposed to recover the property of neutrino number conservation (Bruenn 1985; Liebendorfer et al. 2004). In the 
original version of VERTEX (Rampp & Janka 2002) the problem is overcome by solving both the monochromatic 
moment equations for the neutrino energy density and flux (J, iJ), and the number density and flux {J^ ,%) for 
electron neutrinos and antineutrinos simultaneously, but there are significant drawbacks to this approach: First, the 
mean energy J/ J within a given energy zone [£i-i/2, £^+1/2] is not constrained to the cell-center value Si and may even 
move beyond the cell boundaries if J is treated as an independent quantity. Moreover, the solution of the equations 
by the Newton- Raphson method involves the inversion of square matrix blocks of size 4 x A'g -I- 2 instead of 2 x TV^ 4- 2 
for neutrino energy transport alone, which makes the Newton-Raphson iteration almost eight times more expensive. 

Fortunately, the Doppler terms in the moment equations for neutrino energy transport can be discretized in such 
a way as to guarantee neutrino energy and number conservation simultaneously. For example, in a rather simple 
approach the interface values J2+i/2 be computed iteratively from Ji and Ji-i/2 using the condition 

£i-hl/2 — £i -.n+l -r £i — £j-l/2 ,„+i ,-r,.x 
^«.+1/2^;Y/2 = ^^J^- —^^^-l/2J"-l/2^ (B4) 

SO that the last three terms in Eq. (|B3|) would indeed cancel. However, since this scheme allows large violations of the 
monotonicity condition {Ji < Ji+i/2 < Ji+i or Ji > Ji+1/2 > Ji+i), the resulting algorithm would be rather unstable. 
To construct a stable scheme, we need to impose a weaker condition on the interface fluxes F = ewJ, i.e. we only 
require the source term in the frequency-integrated neutrino number equation to vanish: 

J^tot ~ J^tot\ jn+l £i-l/2 jn+1 rn+l ^£A n mc;^ 
j D " ^ [—^^^+l/2J^+l/2 —W^-l/2J^-l/2 ~ " J = 0" (^5) 

In order to find interface fluxes -Fi+1/2 = £i+i/2Wi+i/2J2+i/2 that fulfil this condition, we first rearrange the sum in 
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Eq. (IB5|) in the following manner: 
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>^t,nt 
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(B6) 
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Next we split the fluxes into "left" and "right" components which are assigned to the adjacent energy zones, 



F- 



i+l/2 — ^ i -r ^ i+i, 

and again rewrite the frequency-integrated Doppler term in the discretized neutrino number equation: 



(B8) 
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(B9) 



If and are chosen such that each of the square brackets in Eq. (|B9I) vanishes, no source term in the frequency- 
integrated neutrino number equation appears, and hence the total lepton number is also conserved^'^. 
The condition we obtain for the "left" and "right" fluxes. 



£i-l 



F'' 



-i+l 



(BIO) 



only fixes the (weighted) sum of and Fl, so we are still free to specify the ratio Fh/Fl to obtain steeper flux 
profiles in the high-energy tail of the spectrum, where Ji decreases rapidly with the zone index i. To this end, we 
parameterize and F^ by a weighting factor 



F^=- 



1 - 

WiAsi 



1 



-J- 



n+l 



(1-6). 



(Bll) 
(B12) 



Here, ^ is defined as a function of the zeroth angular moment j = Jh^c^ /e^ of the neutrino distribution function at 
zone interfaces (obtained as weighted geometric mean of j in the adjacent zones), and a steepness parameter a, 



h-1/2 + •?'i+l/2 



(B13) 



Only the lowest energy zone constitutes an exception, since the condition that F{e = 0) should vanish requires = 1. 
For the steepness parameter, we typically set cr = 1/2 or (7 = 1; both choices result in a second-order accurate scheme^^. 

The rationale for this seemingly complicated procedure for determining the weighting factor can be illustrated as 
follows: First, we specialize to the standard case of a logarithmic grid as used in VERTEX so that Eqs. (|B11[) and 

Strictly speaking, this statement holds only as long as the non-linear moment equations are solved exactly. In practice, the time-critical 
Newton-Raphson iteration is terminated once a specified accuracy in the solution is reached, and a sufficiently low tolerance level must be 
chosen to guarantee long-time conservation of lepton number. 

Second-order accuracy can be proved by means of a Taylor series expansion, but this involves lengthy calculations that are omitted 

here. 

More precisely, the zone interfaces and centers are chosen as follows: £^+1/2 = ^£0-^' (for * G {1, • ■ ■ , -'V}), £-1/2 = Oi and £i = 
5 (ei_i/2 -I- £i+i/2) • Here, A'' is the number of energy zones, Aeg is the width of the first zone, and A > 1 is an adjustable parameter. 
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(|B12I) simplify to, 

Ft^w.,e.,+,/,J]^+'^., (B14) 

i^«=i.,e,_i/2Jr+i(l-e»). (B15) 

In regions where Ji varies mildly with i, ^ is close to 1/2, and hence the total interface flux is approximately given by 
the arithmetic mean of the flux in the adjacent cells, 

1+1/2 - i'i + i'l+i ~ ^ • (Blbj 

On the other hand, in the high-energy tail of the spectrum Ji drops rapidly with increasing z, i.e. J^+i ^ J^. Assuming 
that the spectrum can be approximated by a Fermi-Dirac distribution (with a small chemical potential in the sense 
that jjL <^ £i), Ji is given in terms of the inverse temperature /3 and a normalization factor k as 

J, sa ke-^^'el (B17) 

and the weighting factor is approximately 
In this case the flux components read 



< 1. (B18) 



Ft = w,e,+i/2fce-'5(^'+'^^^-)ef , (B19) 
Ff^^w^e,^i/2ke-^'^el (B20) 

Bearing in mind that varies far more slowly than e^^^ in the steep tail of the spectrum, we obtain 

Ft « «i»e,+i/2fcJ,+i/2, l-F^^i « |u;»£,_i/2fcJ,| « \Ft__^\ , F,+^,2 « F}- (B21) 

for a = 1/2, while for cr = 1 we have, 

Ft ~ WiEi+i/afcJj+i, Ff^ sa WiCj^i/afc Jj -Fi+i/2 ~ {wi + w^+i) £^+1/2^ Jj+i- (B22) 

Thus, in the case of tr = 1/2, the effective interface value of J appearing in the fluxes is given by the weighted 
geometric mean of J in the adjacent zones; for cr = 1 it is twice the minimum value in either zone, which corresponds 
more closely (but not exactly) to the weighted harmonic mean. Since (7 = 1 generally reduces the absolute values of 
the fluxes, the resulting scheme is somewhat more robust than for a = 1/2, i.e. the Newton- Raphson iteration fails to 
convergence less frequently. 

"Harmonic" interpolation with cr = 1 has another welcome side effect in addition to the improved robustness of the 
scheme: The method of geometric interpolation (also in the original scheme of Rampp & Janka, 2002) suffers from 
a mild "kink" instability, that may produce an unphysical dip in the neutrino spectra. In a long-time simulation of 
the neutron star cooling phase with VERTEX-PROMETHEUS, in which our new scheme was switched on about 
1 s after bounce, we found that "harmonic" interpolation eliminates such a dip, as shown in Fig. [151 Because of the 
favourable stability properties and the smoother neutrino spectra, we preferably work with cr = 1, but cr = 1/2 is also 
a viable choice. 
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